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