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: apla88 - Elementsteifigkeitsroutine
40 * ab88 - Berechnung der Matrizen bbi & bsv
41 * 8-Knoten Serendipity Reissner- Mindlin- Platte
42 * 2.1.2010 Rieg
43 ***********************************************************************/
44
45 /***********************************************************************
46 * Fuer UNIX
47 ***********************************************************************/
48 #ifdef FR_UNIX
49 #include <z88r.h>
50 #endif
51
52 /***********************************************************************
53 * Fuer Windows
54 ***********************************************************************/
55 #ifdef FR_WIN
56 #include <z88r.h>
57 #endif
58
59 /***********************************************************************
60 * Fuer Windows und GTK+
61 ***********************************************************************/
62 #ifdef FR_GTKWIN
63 #include <z88r.h>
64 #endif
65
66 /***********************************************************************
67 * Functions
68 ***********************************************************************/
69 int ab88(FR_DOUBLE *det,FR_DOUBLE *r,FR_DOUBLE *s);
70
71 /***********************************************************************
72 * hier beginnt Function apla88
73 ***********************************************************************/
apla88(void)74 int apla88(void)
75 {
76 extern FR_DOUBLEAY se;
77
78 extern FR_DOUBLE xk[],yk[];
79 extern FR_DOUBLE bbi[],bsv[],xx[],dbi[],dsv[],be[],hi[];
80
81 extern FR_DOUBLE emode,rnuee,qparae,riyye;
82
83 extern FR_INT4 intore,ipflag;
84
85 FR_DOUBLE dbbi[10],dbsv[5];
86
87 FR_DOUBLE facbi,facsv,r,s,det,wt,stiffbi,stiffsv,rmok,skf;
88
89 FR_INT4 ne= 24,i,lx,ly,j,k,l;
90
91 int iret;
92
93 /*----------------------------------------------------------------------
94 * Gauss-Legendre Stuetzstellen
95 *---------------------------------------------------------------------*/
96 static FR_DOUBLE xg[17]= { 0.,
97 0., -.5773502691896, -.7745966692415, -.8611363115941,
98 0., +.5773502691896, 0., -.3399810435849,
99 0., 0., +.7745966692415, +.3399810435849,
100 0., 0., 0., +.8611363115941 };
101
102 /*----------------------------------------------------------------------
103 * Gauss-Legendre Integrationsgewichte
104 *---------------------------------------------------------------------*/
105 static FR_DOUBLE wgt[17]= { 0.,
106 2., 1., +.5555555555556, +.3478548451375,
107 0., 1., +.8888888888889, +.6521451548625,
108 0., 0., +.5555555555556, +.6521451548625,
109 0., 0., 0., +.3478548451375 };
110
111 /*----------------------------------------------------------------------
112 * xk und yk umspeichern
113 *---------------------------------------------------------------------*/
114 for(i = 1;i <= 8;i++)
115 {
116 xx[i] = xk[i];
117 xx[8+i]= yk[i];
118 }
119
120 /*----------------------------------------------------------------------
121 * Elastizitaetsmatrix aufstellen: Platten-Biegung
122 *---------------------------------------------------------------------*/
123 facbi = emode*qparae*qparae*qparae/(12.*(1. - rnuee*rnuee));
124 dbi[1]= facbi;
125 dbi[2]= facbi * rnuee;
126 dbi[3]= 0.;
127 dbi[4]= dbi[2];
128 dbi[5]= dbi[1];
129 dbi[6]= 0.;
130 dbi[7]= 0.;
131 dbi[8]= 0.;
132 dbi[9]= facbi * .5 * (1. - rnuee);
133
134 /*----------------------------------------------------------------------
135 * Elastizitaetsmatrix aufstellen: transversale Schubverzerrung
136 *---------------------------------------------------------------------*/
137 rmok= 1;
138 if(ipflag == 1) rmok= 1.; /* Reissner- Mindlin */
139 if(ipflag == 2) rmok= 0.1; /* Schubeinfluss daempfen */
140 if(ipflag == 3) rmok= 0.01; /* Schubeinfluss daempfen */
141 if(ipflag == 4) rmok= 0.001; /* Schubeinfluss daempfen */
142
143 skf= 5./6.; /* Schubkorrekturfaktor */
144
145 facsv= rmok*emode*skf*qparae/(2*(1. + rnuee));
146 dsv[1]= facsv;
147 dsv[2]= 0.;
148 dsv[3]= 0.;
149 dsv[4]= facsv;
150
151 /*----------------------------------------------------------------------
152 * Elementsteifigkeitsmatrix aufstellen
153 *---------------------------------------------------------------------*/
154 for(i = 1;i <= 576;i++)
155 se[i]= 0.;
156
157 for(i = 1;i <= 24;i++)
158 be[i]= 0.;
159
160 for(lx = 1;lx <= intore;lx++)
161 {
162 r= xg[(lx-1)*4 + intore];
163 for(ly = 1;ly <= intore;ly++)
164 {
165 s= xg[(ly-1)*4 + intore];
166
167 /*======================================================================
168 * Matrix b der partiellen Ableitungen & Jacobi Determinante holen
169 *=====================================================================*/
170 iret= ab88(&det,&r,&s);
171 if(iret != 0) return(iret);
172
173 wt= wgt[(lx-1)*4 + intore] * wgt[(ly-1)*4 + intore] * det;
174
175 /*======================================================================
176 * Element- Lastvektor be
177 *=====================================================================*/
178 for(j = 1;j <= 24;j++)
179 {
180 be[j]+= hi[j]*wt*riyye;
181 }
182
183 /*======================================================================
184 * Start Steifigkeitsmatrix
185 *=====================================================================*/
186 for(j = 1;j <= 24;j++)
187 {
188
189 /*======================================================================
190 * Biegeverzerrung: DBBI= B*C fuer Biegung
191 *=====================================================================*/
192 for(k = 1;k <= 3;k++)
193 {
194 dbbi[k]= 0.;
195 for(l = 1;l <= 3;l++)
196 {
197 dbbi[k]+= dbi[(k-1)*3 + l] * bbi[(l-1)*24 + j];
198 }
199 }
200
201 /*======================================================================
202 * Schubverzerrung: DBSV= B*C fuer Schub
203 *=====================================================================*/
204 for(k = 1;k <= 2;k++)
205 {
206 dbsv[k]= 0.;
207 for(l = 1;l <= 2;l++)
208 {
209 dbsv[k]+= dsv[(k-1)*2 + l] * bsv[(l-1)*24 + j];
210 }
211 }
212
213 /*======================================================================
214 * Steifigkeitsmatrix: Die jeweiligen DB's * B und aufsummieren
215 *=====================================================================*/
216 for(i = j;i <= 24;i++)
217 {
218 stiffbi= 0.;
219 stiffsv= 0.;
220
221 for(l = 1;l <= 3;l++)
222 stiffbi+= bbi[(l-1)*24 + i] * dbbi[l];
223
224 for(l = 1;l <= 2;l++)
225 stiffsv+= bsv[(l-1)*24 + i] * dbsv[l];
226
227 se[i+ne*(j-1)]= se[i+ne*(j-1)] + (stiffbi+stiffsv) * wt;
228 }
229 }
230 }
231 }
232
233 /*======================================================================
234 * die andere Haelfte der Steifigkeitsmatrix
235 *=====================================================================*/
236 for(j = 1;j <= 24;j++)
237 {
238 for(i = j;i <= 24;i++)
239 {
240 se[j+ne*(i-1)]= se[i+ne*(j-1)];
241 }
242 }
243
244 return(0);
245 }
246
247 /***********************************************************************
248 * hier beginnt Function ab88
249 ***********************************************************************/
ab88(FR_DOUBLE * det,FR_DOUBLE * r,FR_DOUBLE * s)250 int ab88(FR_DOUBLE *det,FR_DOUBLE *r,FR_DOUBLE *s)
251 {
252 /*---------------------------------------------------------------------
253 * xx geht rein, unveraendert (ex)
254 * bbi und bsv gehen raus, neu (ex)
255 * det geht raus, neu
256 * r,s gehen rein, unveraendert
257 *--------------------------------------------------------------------*/
258
259 extern FR_DOUBLE h[];
260 extern FR_DOUBLE bbi[],bsv[],xx[],p[],hi[];
261
262 FR_DOUBLE xj[5], xji[5]; /* ist 2x2 +1 */
263
264 FR_DOUBLE rp,sp,rm,sm,rqm,sqm,r2,s2,dum;
265
266 FR_INT4 i,j,k,k3;
267
268 /*----------------------------------------------------------------------
269 * Klammern der Formfunktionen belegen
270 *---------------------------------------------------------------------*/
271 rp= 1. + (*r);
272 sp= 1. + (*s);
273 rm= 1. - (*r);
274 sm= 1. - (*s);
275 rqm= 1. - (*r)*(*r);
276 sqm= 1. - (*s)*(*s);
277 r2= 2. * (*r);
278 s2= 2. * (*s);
279
280 /*----------------------------------------------------------------------
281 * Formfunktionen
282 *---------------------------------------------------------------------*/
283 h[1]= .25 *(rp*sp - rqm*sp - sqm*rp);
284 h[2]= .25 *(rm*sp - rqm*sp - sqm*rm);
285 h[3]= .25 *(rm*sm - sqm*rm - rqm*sm);
286 h[4]= .25 *(rp*sm - rqm*sm - sqm*rp);
287 h[5]= .5 *rqm*sp;
288 h[6]= .5 *sqm*rm;
289 h[7]= .5 *rqm*sm;
290 h[8]= .5 *sqm*rp;
291
292 /*----------------------------------------------------------------------
293 * Partielle Ableitung der Formfunktionen nach r
294 *---------------------------------------------------------------------*/
295 p[1]= .25 *(sp + r2*sp -sqm);
296 p[2]= .25 *((-sp) + r2*sp + sqm);
297 p[3]= .25 *((-sm) + sqm + r2*sm);
298 p[4]= .25 *(sm + r2*sm - sqm);
299 p[5]= .5 *(-r2)*sp;
300 p[6]= (-.5 )*sqm;
301 p[7]= .5 *(-r2)*sm;
302 p[8]= .5 *sqm;
303
304 /*----------------------------------------------------------------------
305 * Partielle Ableitung der Formfunktionen nach s
306 *---------------------------------------------------------------------*/
307 p[9] = .25 *(rp - rqm + s2*rp);
308 p[10]= .25 *(rm - rqm + s2*rm);
309 p[11]= .25 *((-rm) + s2*rm + rqm);
310 p[12]= .25 *((-rp) + rqm + s2*rp);
311 p[13]= .5 *rqm;
312 p[14]= .5 *(-s2)*rm;
313 p[15]= (-.5 )*rqm;
314 p[16]= .5 *(-s2)*rp;
315
316 /*----------------------------------------------------------------------
317 * Jacobi-Matrix am Punkt (r,s) entwickeln
318 *---------------------------------------------------------------------*/
319 for(i = 1;i <= 2;i++)
320 {
321 for(j = 1;j <= 2;j++)
322 {
323 dum= 0.;
324 for(k = 1;k <= 8;k++)
325 {
326 dum+= p[(i-1)*8 + k] * xx[(j-1)*8 + k];
327 }
328 xj[(i-1)*2 + j]= dum;
329 }
330 }
331
332 /*----------------------------------------------------------------------
333 * Jacobi-Determinante am Punkt (r,s) entwickeln
334 *---------------------------------------------------------------------*/
335 (*det)= xj[1] * xj[4] - xj[3] * xj[2];
336
337 if((*det) < 0.00000001)
338 return(AL_JACNEG);
339
340 /*----------------------------------------------------------------------
341 * Berechnung der inversen Jacobi-Matrix
342 *---------------------------------------------------------------------*/
343 dum= 1./(*det);
344
345 xji[1]= xj[4] * dum;
346 xji[2]= (-xj[2]) * dum;
347 xji[3]= (-xj[3]) * dum;
348 xji[4]= xj[1] * dum;
349
350 /*----------------------------------------------------------------------
351 * Entwickeln der Matrix bbi fuer Biegung
352 *---------------------------------------------------------------------*/
353 for(i = 1;i <= 3*24;i++)
354 bbi[i]= 0.;
355
356 k3= 0;
357
358 for(k = 1;k <= 8;k++)
359 {
360 k3+= 3;
361
362 for(i = 1;i <= 2;i++)
363 {
364 bbi[ k3 ]+= xji[ i] * p[(i-1)*8+k];
365 bbi[24+k3-1]-= xji[2+i] * p[(i-1)*8+k];
366 }
367 bbi[48+k3 ]= -bbi[24+k3-1];
368 bbi[48+k3-1]= -bbi[ k3 ];
369 }
370
371 /*----------------------------------------------------------------------
372 * Entwickeln der Matrix bsv fuer Schub
373 *---------------------------------------------------------------------*/
374 for(i = 1;i <= 2*24;i++)
375 bsv[i]= 0.;
376
377 k3= 0;
378
379 for(k = 1;k <= 8;k++)
380 {
381 k3+= 3;
382
383 for(i = 1;i <= 2;i++)
384 {
385 bsv[ k3-2]+= xji[2+i] * p[(i-1)*8+k];
386 bsv[24+k3-2]+= xji[ i] * p[(i-1)*8+k];
387 }
388 bsv[ k3-1]= -h[k];
389 bsv[24+k3 ]= h[k];
390 }
391
392 /*----------------------------------------------------------------------
393 * Entwickeln der Formfunktionen fuer den Lastvektor be
394 *---------------------------------------------------------------------*/
395 for(i = 1;i <= 24;i++)
396 hi[i]= 0.;
397
398 k3= 1;
399
400 for(k = 1;k <= 8;k++)
401 {
402 hi[k3]= h[k];
403 k3+= 3;
404 }
405
406 return(0);
407 }
408