1 /***********************************************************************
2 *
3 * ***** *** ***
4 * * * * * *
5 * * *** ***
6 * * * * * *
7 * ***** *** ***
8 *
9 * A FREE Finite Elements Analysis Program in ANSI C for the Windows &
10 * the 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: sshf88 - Elementsteifigkeitsroutine
37 * ruft: ibs88 - Berechnung Matrix b - Schale
38 * ruft: sss88 - Berechnung Matrix b - Platte
39 * 24.5.2012 Rieg
40 ***********************************************************************/
41
42 /***********************************************************************
43 * Fuer UNIX
44 ***********************************************************************/
45 #ifdef FR_UNIX
46 #include <z88r.h>
47 #include <stdio.h> /* fprintf */
48 #endif
49
50 /***********************************************************************
51 * Fuer Windows
52 ***********************************************************************/
53 #ifdef FR_WIN
54 #include <z88r.h>
55 #include <stdio.h> /* fprintf */
56 #endif
57
58 /***********************************************************************
59 * Fuer Windows und GTK+
60 ***********************************************************************/
61 #ifdef FR_GTKWIN
62 #include <z88r.h>
63 #include <stdio.h> /* fprintf */
64 #endif
65
66 /***********************************************************************
67 * Formate
68 ***********************************************************************/
69 #define NL "\n"
70
71 #ifdef FR_XDOUB
72 #define P10E " %+#10.2lE"
73 #define P11E " %+#11.3lE"
74 #endif
75
76 #ifdef FR_XQUAD
77 #define P10E " %+#10.2LE"
78 #define P11E " %+#11.3LE"
79 #endif
80
81 /***********************************************************************
82 * Functions
83 ***********************************************************************/
84 int ibs88(FR_DOUBLE *det,FR_DOUBLE *r,FR_DOUBLE *s,
85 FR_DOUBLE *xbar,FR_INT4 *ktyp);
86 int sss88(FR_DOUBLE *det,FR_DOUBLE *r,FR_DOUBLE *s);
87 FR_DOUBLE platgh(FR_DOUBLE sig[]);
88 FR_DOUBLE platnh(FR_DOUBLE sig[]);
89 FR_DOUBLE platsh(FR_DOUBLE sig[]);
90
91 /***********************************************************************
92 * hier beginnt Function sshf88
93 ***********************************************************************/
sshf88(void)94 int sshf88(void)
95 {
96 extern FILE *fo3,*fo5;
97
98 extern FR_DOUBLEAY smw;
99 extern FR_DOUBLEAY gmw;
100 extern FR_DOUBLEAY sigvku;
101 extern FR_DOUBLEAY tm;
102 extern FR_DOUBLEAY tmt;
103
104 extern FR_INT4AY jsm;
105 extern FR_INT4AY koi;
106 extern FR_INT4AY koffs;
107
108 extern FR_DOUBLE ul[];
109 extern FR_DOUBLE h[];
110 extern FR_DOUBLE xk[],yk[],zk[];
111 extern FR_DOUBLE b[],xx[],ds[],dp[];
112 extern FR_DOUBLE xc[],yc[],zc[];
113
114 extern FR_DOUBLE emode,rnuee,qparae;
115
116 extern FR_INT4 ktyp,ninto,kflag,isflag,kc,jpri,ifnili,ngau,ihflag;
117
118 FR_DOUBLE ulc[37],uls[13],ulp[19];
119 FR_DOUBLE a21x,a21y,a21z,rl21,a31x,a31y,a31z;
120 FR_DOUBLE azsx,azsy,azsz,rlzs,aysx,aysy,aysz,rlys;
121 FR_DOUBLE vxx,vxy,vxz,vyx,vyy,vyz,vzx,vzy,vzz;
122
123 FR_DOUBLE sigs[6],sigp[6],rsig[6];
124 FR_DOUBLE eps[6],sig[7],rv[4];
125
126 FR_DOUBLE facesz,facasz,facbi,facsv,rmok,skf,fque,fmom;
127 FR_DOUBLE r,s,xs,ys,zs,xsc,ysc,zsc,det,xbar,sum;
128 FR_DOUBLE sigv,ax,ay,az,rmin;
129
130 FR_INT4 jp[4];
131
132 FR_INT4 i,lx,j,k,jk,inc,jnc;
133
134 int iret;
135
136 /*----------------------------------------------------------------------
137 * Gauss-Legendre Stuetzstellen fuer r
138 *---------------------------------------------------------------------*/
139 static FR_DOUBLE rg[40]=
140 {
141 0.,0.,0.,0.,0.,0.,0., /* Elemente 0 - 6 leer */
142 0.1666666666667, /* intore = 3, 1.Ele Start bei i=7 */
143 0.6666666666667,
144 0.1666666666667,
145 0.,0.,0.,0.,0., /* Elemente 10-14 leer */
146 0.1012865073235, /* intore = 7, 1.Ele Start bei i=15 */
147 0.7974269853531,
148 0.1012865073235,
149 0.4701420641051,
150 0.4701420641051,
151 0.0597158717898,
152 0.3333333333333,
153 0.,0.,0.,0.,0., /* Elemente 22-26 leer */
154 0.0651301029022, /* intore =13, 1.Ele Start bei i=27 */
155 0.8697397941956,
156 0.0651301029022,
157 0.3128654960049,
158 0.6384441885698,
159 0.0486903154253,
160 0.6384441885698,
161 0.3128654960049,
162 0.0486903154253,
163 0.2603459660790,
164 0.4793080678419,
165 0.2603459660790,
166 0.3333333333333
167 };
168
169 /*----------------------------------------------------------------------
170 * Gauss-Legendre Stuetzstellen fuer s
171 *---------------------------------------------------------------------*/
172 static FR_DOUBLE sg[40]=
173 {
174 0.,0.,0.,0.,0.,0.,0., /* Elemente 0 - 6 leer */
175 0.1666666666667, /* intore = 3, 1.Ele Start bei i=7 */
176 0.1666666666667,
177 0.6666666666667,
178 0.,0.,0.,0.,0., /* Elemente 10-14 leer */
179 0.1012865073235, /* intore = 7, 1.Ele Start bei i=15 */
180 0.1012865073235,
181 0.7974269853531,
182 0.0597158717898,
183 0.4701420641051,
184 0.4701420641051,
185 0.3333333333333,
186 0.,0.,0.,0.,0., /* Elemente 22-26 leer */
187 0.0651301029022, /* intore =13, 1.Ele Start bei i=27 */
188 0.0651301029022,
189 0.8697397941956,
190 0.0486903154253,
191 0.3128654960049,
192 0.6384441885698,
193 0.0486903154253,
194 0.6384441885698,
195 0.3128654960049,
196 0.2603459660790,
197 0.2603459660790,
198 0.4793080678419,
199 0.3333333333333
200 };
201
202 /*----------------------------------------------------------------------
203 * Gauss-Legendre Stuetzstellen, fix fuer 3 Punkte
204 *---------------------------------------------------------------------*/
205 static FR_DOUBLE xgo[4]= { 0.,
206 0.1666666666667,
207 0.6666666666667,
208 0.1666666666667};
209
210 static FR_DOUBLE ygo[4]= { 0.,
211 0.1666666666667,
212 0.1666666666667,
213 0.6666666666667};
214
215 /*----------------------------------------------------------------------
216 * Natuerliche Koordinaten der Eckknoten
217 *---------------------------------------------------------------------*/
218 static FR_DOUBLE rkr[4]= { 0.,
219 0., 1., 0. };
220 static FR_DOUBLE rks[4]= { 0.,
221 0., 0., 1. };
222
223 /*----------------------------------------------------------------------
224 * globale Koordinaten in lokale Koordinaten umrechnen
225 *---------------------------------------------------------------------*/
226 /*======================================================================
227 * Vektor 2-1 spannt die lokale x-Achse x' auf
228 *=====================================================================*/
229 a21x= xk[2] - xk[1];
230 a21y= yk[2] - yk[1];
231 a21z= zk[2] - zk[1];
232 rl21 = sqrt(a21x*a21x + a21y*a21y + a21z*a21z);
233
234 /*======================================================================
235 * Vektor 3-1 wird fuer Kreuzprodukt azs gebraucht
236 *=====================================================================*/
237 a31x= xk[3]-xk[1];
238 a31y= yk[3]-yk[1];
239 a31z= zk[3]-zk[1];
240
241 /*======================================================================
242 * azs ist Kreuzprodukt aus a21 und a31 und bildet die lokale Achse z'
243 *=====================================================================*/
244 azsx= a21y*a31z - a21z*a31y;
245 azsy= a21z*a31x - a21x*a31z;
246 azsz= a21x*a31y - a21y*a31x;
247 rlzs = sqrt(azsx*azsx + azsy*azsy + azsz*azsz);
248
249 /*======================================================================
250 * ays ist Kreuzprodukt aus azs und a21 und bildet die lokale Achse y'
251 *=====================================================================*/
252 aysx= azsy*a21z - azsz*a21y;
253 aysy= azsz*a21x - azsx*a21z;
254 aysz= azsx*a21y - azsy*a21x;
255 rlys = sqrt(aysx*aysx + aysy*aysy + aysz*aysz);
256
257 /*======================================================================
258 * Richtungscosinus x'x, x'y und x'z
259 *=====================================================================*/
260 vxx=a21x/rl21;
261 vxy=a21y/rl21;
262 vxz=a21z/rl21;
263
264 /*======================================================================
265 * Richtungscosinus y'x, y'y und y'z
266 *=====================================================================*/
267 vyx=aysx/rlys;
268 vyy=aysy/rlys;
269 vyz=aysz/rlys;
270
271 /*======================================================================
272 * Richtungscosinus z'x, z'y und z'z
273 *=====================================================================*/
274 vzx=azsx/rlzs;
275 vzy=azsy/rlzs;
276 vzz=azsz/rlzs;
277
278 /*======================================================================
279 * xk auf xc zwischensichern, damit xk lokal werden kann
280 *=====================================================================*/
281 for(i= 1; i <= 6; i++)
282 {
283 xc[i]= xk[i];
284 yc[i]= yk[i];
285 zc[i]= zk[i];
286 }
287
288 /*======================================================================
289 * nun Koordinaten transformieren, xk,yk,zk sind nun lokal
290 *=====================================================================*/
291 for(i= 1; i <= 6; i++)
292 {
293 xk[i]= vxx*xc[i] + vxy*yc[i] + vxz*zc[i];
294 yk[i]= vyx*xc[i] + vyy*yc[i] + vyz*zc[i];
295 zk[i]= vzx*xc[i] + vzy*yc[i] + vzz*zc[i];
296 }
297
298 /*----------------------------------------------------------------------
299 * xk und yk auf xx umspeichern, zk fuer die spaetere Ruecktransforma.
300 *---------------------------------------------------------------------*/
301 for(i = 1;i <= 6;i++)
302 {
303 xx[ i]= xk[i];
304 xx[ 6+i]= yk[i];
305 xx[12+i]= zk[i];
306 }
307
308 /*----------------------------------------------------------------------
309 * die Transformationsmatrix tm und die Transponierte tmt aufstellen
310 *---------------------------------------------------------------------*/
311 /*======================================================================
312 * Transformationsmatrix tm
313 *=====================================================================*/
314 for(i = 1; i <= 1296; i++)
315 tm[i]= 0.;
316
317 for(i = 1; i <= 1296; i++)
318 tmt[i]= 0.;
319
320 inc= 0;
321
322 for(i = 1; i <= 12; i++)
323 {
324 tm[inc*36+inc +1]= vxx;
325 tm[inc*36+inc +2]= vxy;
326 tm[inc*36+inc +3]= vxz;
327 tm[inc*36+inc+36+1]= vyx;
328 tm[inc*36+inc+36+2]= vyy;
329 tm[inc*36+inc+36+3]= vyz;
330 tm[inc*36+inc+72+1]= vzx;
331 tm[inc*36+inc+72+2]= vzy;
332 tm[inc*36+inc+72+3]= vzz;
333 inc+= 3;
334 }
335
336 /*======================================================================
337 * Transponierte tmt der Transformationsmatrix tm
338 *=====================================================================*/
339 for(i= 1; i <= 36; i++)
340 {
341 for(j= 1; j <= 36; j++)
342 {
343 tmt[(i-1)*36+j]= tm[(j-1)*36+i];
344 }
345 }
346
347 /*----------------------------------------------------------------------
348 * die globalen Verschiebungen zwischensichern
349 *---------------------------------------------------------------------*/
350 for(i = 1; i <= 36; i++)
351 ulc[i]= ul[i];
352
353 /*----------------------------------------------------------------------
354 * Matrix-Vektorprodukt: tm * ulc(global) = ul(lokal)
355 *---------------------------------------------------------------------*/
356 for(i = 1;i <= 36;i++)
357 {
358 sum= 0.;
359 for(k = 1; k <= 36; k++)
360 {
361 sum+= tm[(i-1)*36+k] * ulc[k];
362 }
363 ul[i]= sum; /* ul ist jetzt lokal */
364 }
365
366 /*----------------------------------------------------------------------
367 * der Scheibenanteil von ul(lokal)
368 *---------------------------------------------------------------------*/
369 inc= 1;
370 jnc= 1;
371 for(i = 1;i <= 6;i++)
372 {
373 uls[jnc ]= ul[inc ];
374 uls[jnc+1]= ul[inc+1];
375 jnc+=2;
376 inc+=6;
377 }
378
379 /*----------------------------------------------------------------------
380 * der Plattenanteil von ul(lokal)
381 *---------------------------------------------------------------------*/
382 inc= 3;
383 jnc= 1;
384 for(i = 1;i <= 6;i++)
385 {
386 ulp[jnc ]= ul[inc ];
387 ulp[jnc+1]= ul[inc+1];
388 ulp[jnc+2]= ul[inc+2];
389 jnc+=3;
390 inc+=6;
391 }
392
393 /*----------------------------------------------------------------------
394 * Scheibe: Materialmatrix aufstellen
395 *---------------------------------------------------------------------*/
396 for(i = 1;i <= 36;i++)
397 ds[i]= 0.;
398
399 ktyp= 2;
400
401 facesz= emode/(1. - rnuee*rnuee);
402 facasz= emode*(1. - rnuee)/( (1. + rnuee)*(1. - 2*rnuee) );
403
404 ds[1] = facesz;
405 ds[5] = facesz * rnuee;
406 ds[9] = 0.;
407 ds[2] = ds[5];
408 ds[6] = facesz;
409 ds[10]= 0.;
410 ds[3] = 0.;
411 ds[7] = 0.;
412 ds[11]= facesz * .5 * (1. - rnuee);
413
414 /*----------------------------------------------------------------------
415 * Platte: Materialmatrix aufstellen
416 *---------------------------------------------------------------------*/
417 for(i = 1;i <= 36;i++)
418 dp[i]= 0.;
419
420 facbi = emode*qparae*qparae*qparae/(12.*(1. - rnuee*rnuee));
421
422 dp[1] = facbi;
423 dp[2] = facbi * rnuee;
424
425 dp[6] = dp[2];
426 dp[7] = dp[1];
427
428 dp[13]= facbi * .5 * (1. - rnuee);
429
430 if (ihflag == 1) rmok= 1.; /* Reissner- Mindlin */
431 else if(ihflag == 2) rmok= 0.01; /* Schubeinfluss daempfen */
432 else rmok= 1.; /* Default */
433
434 skf= 5./6.; /* Schubkorrekturfaktor */
435
436 facsv= rmok*emode*skf*qparae/(2*(1. + rnuee));
437
438 dp[19]= facsv;
439 dp[25]= facsv;
440
441 /*----------------------------------------------------------------------
442 * Spannungen in den Gauss-Punkten berechnen
443 *---------------------------------------------------------------------*/
444 if(ninto > 0)
445 {
446
447 /*======================================================================
448 * ninto > 0: Spannungen in den Gauss-Punkten berechnen, variabel
449 *=====================================================================*/
450 for(lx = 1;lx <= ninto;lx++)
451 {
452 r= rg[lx+2*ninto];
453 s= sg[lx+2*ninto];
454
455 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
456 * ninto > 0: Scheibe: Matrix b & Formfunktionen holen
457 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
458 iret= ibs88(&det,&r,&s,&xbar,&ktyp);
459 if(iret != 0) return(iret);
460
461 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
462 * ninto > 0: Scheibe: Verzerrungen berechnen
463 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
464 for(k = 1;k <= 3;k++)
465 {
466 eps[k]= 0.;
467 for(j = 1;j <= 12;j++)
468 {
469 eps[k]+= b[(k-1)*12 + j] * uls[j];
470 }
471 }
472
473 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
474 * ninto > 0: Scheibe: Spannungen berechnen
475 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
476 for(k = 1;k <= 3;k++)
477 {
478 sigs[k]= 0.;
479 for(j = 1;j <= 3;j++)
480 {
481 sigs[k]+= ds[(k-1)*4 + j] * eps[j];
482 }
483 }
484
485 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
486 * ninto > 0: Platte: Matrix b & Formfunktionen holen
487 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
488 iret= sss88(&det,&r,&s);
489 if(iret != 0) return(iret);
490
491 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
492 * ninto > 0: Platte: Verzerrungen berechnen
493 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
494 for(k = 1;k <= 5;k++)
495 {
496 eps[k]= 0.;
497 for(j = 1;j <= 18;j++)
498 {
499 eps[k]+= b[(k-1)*18 + j] * ulp[j];
500 }
501 }
502
503 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
504 * ninto > 0: Platte: Spannungen berechnen: Biegemomente und Querkraefte
505 * pro Laenge
506 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
507 for(k = 1;k <= 5;k++)
508 {
509 sigp[k]= 0.;
510 for(j = 1;j <= 5;j++)
511 {
512 sigp[k]+= dp[(k-1)*5 + j] * eps[j];
513 }
514 }
515
516 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
517 * ninto > 0: Platte: echte Spannungen berechnen
518 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
519 fmom= 12.*(qparae/2.)/(qparae*qparae*qparae);
520
521 for(k = 1;k <= 3;k++)
522 rsig[k]= sigp[k]*fmom;
523
524 fque= 3./2. /qparae;
525
526 for(k = 4;k <= 5;k++)
527 rsig[k]= sigp[k]*fque;
528
529 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
530 * ninto > 0: Schale: aus Scheibe und Platte die Spannungen zusammenbauen
531 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
532 sig[1]= sigs[1]+rsig[1]; /* SIGXX */
533 sig[2]= sigs[2]+rsig[2]; /* SIGYY */
534 sig[3]= sigs[3]+rsig[3]; /* TAUXY */
535 sig[4]= 0; /* TAUXZ,d.h. rsig[4],vernachlaessigen */
536 sig[5]= 0; /* TAUYZ,d.h. rsig[5],vernachlaessigen */
537 sig[6]= 0.;
538
539 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
540 * ninto > 0: Schale: Integrationspunkte in echte Koordinaten umrechnen
541 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
542 xsc= 0.;
543 ysc= 0.;
544 zsc= 0.;
545
546 /* von r,t auf lokal */
547 for(j = 1;j <= 6;j++)
548 {
549 xsc+= h[j] * xx[ j];
550 ysc+= h[j] * xx[6 +j];
551 zsc+= h[j] * xx[12+j];
552 }
553
554 /* von lokal auf global */
555 xs= tmt[ 1]*xsc + tmt[ 2]*ysc + tmt[ 3]*zsc;
556 ys= tmt[37]*xsc + tmt[38]*ysc + tmt[39]*zsc;
557 zs= tmt[73]*xsc + tmt[74]*ysc + tmt[75]*zsc;
558
559 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
560 * ninto > 0: Schale: Spannungen ausschreiben
561 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
562 if(isflag == 0) /* keine VglSpa.*/
563 {
564 if(ifnili == 0)
565 {
566 fprintf(fo3,NL P10E P10E P10E P11E P11E P11E
567 ,xs,ys,zs,sig[1],sig[2],sig[3]);
568 }
569 } /* e if ohne VglS */
570
571 if(isflag == 1 || isflag == 2 || isflag == 3) /* GEH,NH,SH */
572 {
573 if(isflag == 1) sigv= platgh(sig);
574 if(isflag == 2) sigv= platnh(sig);
575 if(isflag == 3) sigv= platsh(sig);
576
577 ngau++;
578 sigvku[ngau]+= sigv;
579
580 if(ifnili == 0)
581 {
582 fprintf(fo3,NL P10E P10E P10E P11E P11E P11E P10E
583 ,xs,ys,zs,sig[1],sig[2],sig[3],sigv);
584 }
585
586 if(jpri == 1)
587 {
588 fprintf(fo5,NL P11E P11E P11E P11E,xs,ys,zs,sigvku[ngau]);
589 }
590 gmw[kc]+= sigv;
591 } /* e if mit GEH */
592 } /* e for */
593 gmw[kc]/= ninto; /* Mittelwert berechnen */
594
595 /*======================================================================
596 * ninto > 0: die Eckpunkte berechnen
597 *=====================================================================*/
598 for(lx = 1;lx <= 3;lx++)
599 {
600 r= rkr[lx];
601 s= rks[lx];
602
603 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
604 * ninto > 0: Matrix b der partiellen Ableitungen & Formfunktionen holen:
605 * egal, ob * Scheiben- oder Plattenformfunktionen ibs88 bzw. sss88
606 * (identisch, * b nicht!), weil nur h[] gebraucht wird.
607 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
608 iret= ibs88(&det,&r,&s,&xbar,&ktyp);
609 if(iret != 0) return(iret);
610
611 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
612 * ninto > 0: Schale: Eckpunkte in echte Koordinaten umrechnen
613 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
614 xsc= 0.;
615 ysc= 0.;
616 zsc= 0.;
617
618 /* von r,t auf lokal */
619 for(j = 1;j <= 6;j++)
620 {
621 xsc+= h[j] * xx[ j];
622 ysc+= h[j] * xx[6 +j];
623 zsc+= h[j] * xx[12+j];
624 }
625
626 /* von lokal auf global */
627 xs= tmt[ 1]*xsc + tmt[ 2]*ysc + tmt[ 3]*zsc;
628 ys= tmt[37]*xsc + tmt[38]*ysc + tmt[39]*zsc;
629 zs= tmt[73]*xsc + tmt[74]*ysc + tmt[75]*zsc;
630
631 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
632 * ninto > 0: welcher Knoten ist's wirklich?
633 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
634 for(j = 1;j <= 3;j++)
635 {
636 ax = xc[j] - xs;
637 ay = yc[j] - ys;
638 az = zc[j] - zs;
639 rv[j]= FR_SQRT(ax*ax + ay*ay + az*az);
640 }
641
642 rmin= 1e88;
643 for(j = 1;j <= 3;j++)
644 {
645 if(rv[j] < rmin)
646 {
647 rmin= rv[j];
648 jk= j;
649 }
650 }
651
652 jp[lx]= jk;
653 }
654
655 /*======================================================================
656 * ninto > 0: Spannungen in den Gauss-Punkten berechnen, fix fuer Z88O
657 *=====================================================================*/
658 for(lx = 1;lx <= 3;lx++)
659 {
660 r= xgo[lx];
661 s= ygo[lx];
662
663 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
664 * ninto > 0: Scheibe: Matrix b & Formfunktionen holen
665 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
666 iret= ibs88(&det,&r,&s,&xbar,&ktyp);
667 if(iret != 0) return(iret);
668
669 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
670 * ninto > 0: Scheibe: Verzerrungen berechnen
671 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
672 for(k = 1;k <= 3;k++)
673 {
674 eps[k]= 0.;
675 for(j = 1;j <= 12;j++)
676 {
677 eps[k]+= b[(k-1)*12 + j] * uls[j];
678 }
679 }
680
681 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
682 * ninto > 0: Scheibe: Spannungen berechnen
683 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
684 for(k = 1;k <= 3;k++)
685 {
686 sigs[k]= 0.;
687 for(j = 1;j <= 3;j++)
688 {
689 sigs[k]+= ds[(k-1)*4 + j] * eps[j];
690 }
691 }
692
693 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
694 * ninto > 0: Platte: Matrix b & Formfunktionen holen
695 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
696 iret= sss88(&det,&r,&s);
697 if(iret != 0) return(iret);
698
699 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
700 * ninto > 0: Platte: Verzerrungen berechnen
701 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
702 for(k = 1;k <= 5;k++)
703 {
704 eps[k]= 0.;
705 for(j = 1;j <= 18;j++)
706 {
707 eps[k]+= b[(k-1)*18 + j] * ulp[j];
708 }
709 }
710
711 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
712 * ninto > 0: Platte: Spannungen berechnen: Biegemomente und Querkraefte
713 * pro Laenge
714 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
715 for(k = 1;k <= 5;k++)
716 {
717 sigp[k]= 0.;
718 for(j = 1;j <= 5;j++)
719 {
720 sigp[k]+= dp[(k-1)*5 + j] * eps[j];
721 }
722 }
723
724 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
725 * ninto > 0: Platte: echte Spannungen berechnen
726 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
727 fmom= 12.*(qparae/2.)/(qparae*qparae*qparae);
728
729 for(k = 1;k <= 3;k++)
730 rsig[k]= sigp[k]*fmom;
731
732 fque= 3./2. /qparae;
733
734 for(k = 4;k <= 5;k++)
735 rsig[k]= sigp[k]*fque;
736
737 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
738 * ninto > 0: Schale: aus Scheibe und Platte die Spannungen zusammenbauen
739 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
740 sig[1]= sigs[1]+rsig[1]; /* SIGXX */
741 sig[2]= sigs[2]+rsig[2]; /* SIGYY */
742 sig[3]= sigs[3]+rsig[3]; /* TAUXY */
743 sig[4]= 0; /* TAUXZ, ansich rsig[4] */
744 sig[5]= 0; /* TAUYZ, ansich rsig[5] */
745 sig[6]= 0.;
746
747 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
748 * ninto > 0: Schale: Vergleichsspannungen aufaddieren
749 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
750 if(isflag == 1) sigv= platgh(sig);
751 if(isflag == 2) sigv= platnh(sig);
752 if(isflag == 3) sigv= platsh(sig);
753
754 smw[koi[koffs[kc]+jp[lx]-1]]+= sigv;
755 jsm[koi[koffs[kc]+jp[lx]-1]]++;
756 } /* Ende Z88O */
757 } /* Ende Gausspunkte variabel, d.h. ninto > 0 */
758
759 /*----------------------------------------------------------------------
760 * ninto = 0: Spannungen in den Eckknoten berechnen
761 *---------------------------------------------------------------------*/
762 if(ninto == 0)
763 {
764 for(lx = 1;lx <= 3;lx++)
765 {
766 r= rkr[lx];
767 s= rks[lx];
768
769 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
770 * ninto = 0: Scheibe: Matrix b & Formfunktionen holen
771 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
772 iret= ibs88(&det,&r,&s,&xbar,&ktyp);
773 if(iret != 0) return(iret);
774
775 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
776 * ninto = 0: Scheibe: Verzerrungen berechnen
777 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
778 for(k = 1;k <= 3;k++)
779 {
780 eps[k]= 0.;
781 for(j = 1;j <= 12;j++)
782 {
783 eps[k]+= b[(k-1)*12 + j] * uls[j];
784 }
785 }
786
787 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
788 * ninto = 0: Scheibe: Spannungen berechnen
789 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
790 for(k = 1;k <= 3;k++)
791 {
792 sigs[k]= 0.;
793 for(j = 1;j <= 3;j++)
794 {
795 sigs[k]+= ds[(k-1)*4 + j] * eps[j];
796 }
797 }
798
799 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
800 * ninto = 0: Platte: Matrix b & Formfunktionen holen
801 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
802 iret= sss88(&det,&r,&s);
803 if(iret != 0) return(iret);
804
805 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
806 * ninto = 0: Platte: Verzerrungen berechnen
807 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
808 for(k = 1;k <= 5;k++)
809 {
810 eps[k]= 0.;
811 for(j = 1;j <= 18;j++)
812 {
813 eps[k]+= b[(k-1)*18 + j] * ulp[j];
814 }
815 }
816
817 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
818 * ninto = 0: Platte: Spannungen berechnen: Biegemomente und Querkraefte
819 * pro Laenge
820 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
821 for(k = 1;k <= 5;k++)
822 {
823 sigp[k]= 0.;
824 for(j = 1;j <= 5;j++)
825 {
826 sigp[k]+= dp[(k-1)*5 + j] * eps[j];
827 }
828 }
829
830 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
831 * Platte: echte Spannungen berechnen
832 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
833 fmom= 12.*(qparae/2.)/(qparae*qparae*qparae);
834
835 for(k = 1;k <= 3;k++)
836 rsig[k]= sigp[k]*fmom;
837
838 fque= 3./2. /qparae;
839
840 for(k = 4;k <= 5;k++)
841 rsig[k]= sigp[k]*fque;
842
843 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
844 * aus Scheibe und Platte die Spannungen zusammenbauen
845 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
846 sig[1]= sigs[1]+rsig[1]; /* SIGXX */
847 sig[2]= sigs[2]+rsig[2]; /* SIGYY */
848 sig[3]= sigs[3]+rsig[3]; /* TAUXY */
849 sig[4]= 0; /* TAUXZ,d.h. rsig[4],vernachlaessigen */
850 sig[5]= 0; /* TAUYZ,d.h. rsig[5],vernachlaessigen */
851 sig[6]= 0.;
852
853 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
854 * ninto = 0: Schale: Eckpunkte in echte Koordinaten umrechnen
855 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
856 xsc= 0.;
857 ysc= 0.;
858 zsc= 0.;
859
860 /* von r,t auf lokal */
861 for(j = 1;j <= 6;j++)
862 {
863 xsc+= h[j] * xx[ j];
864 ysc+= h[j] * xx[6 +j];
865 zsc+= h[j] * xx[12+j];
866 }
867
868 /* von lokal auf global */
869 xs= tmt[ 1]*xsc + tmt[ 2]*ysc + tmt[ 3]*zsc;
870 ys= tmt[37]*xsc + tmt[38]*ysc + tmt[39]*zsc;
871 zs= tmt[73]*xsc + tmt[74]*ysc + tmt[75]*zsc;
872
873 /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
874 * ninto = 0: Schale: Spannungen ausschreiben
875 *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
876 if(ifnili == 0)
877 {
878 fprintf(fo3,NL P10E P10E P10E P11E P11E P11E
879 ,xs,ys,zs,sig[1],sig[2],sig[3]);
880 }
881 }
882 } /* Ende ninto = 0 */
883
884 return(0);
885 }
886
887
888