/*********************************************************************** * * ***** *** *** * * * * * * * * *** *** * * * * * * * ***** *** *** * * A FREE Finite Elements Analysis Program in ANSI C for the UNIX OS. * * Composed and edited and copyright by * Professor Dr.-Ing. Frank Rieg, University of Bayreuth, Germany * * eMail: * frank.rieg@uni-bayreuth.de * dr.frank.rieg@t-online.de * * V10.0 December 12, 2001 * * Z88 should compile and run under any UNIX OS and Motif 2.0. * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2, or (at your option) * any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; see the file COPYING. If not, write to * the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. ***********************************************************************/ /*********************************************************************** * z88e.c * 2.9.2002 Rieg ***********************************************************************/ /*********************************************************************** * Fuer UNIX ***********************************************************************/ #ifdef FR_UNIX #include #include /* FILE,fclose */ #endif /*********************************************************************** * Functions ***********************************************************************/ int who88e(void); int dyn88e(void); int wrim88e(FR_INT4,int); int wlog88e(FR_INT4,int); int wtyp88e(FR_INT4,FR_INT4); int ale88e(int); int riy88(void); int lqua88(void); int balk88(void); int shei88(void); int stab88(void); int well88(void); int ring88(void); int qshe88(void); int esta88(void); int hexa88(void); int cshe88(void); int ebal88(void); int isod88(void); int tetr88(void); int spur88(void); int spla88(void); int hpla88(void); int apla88(void); void stop88e(void); int forc88(void); int lan88e(void); /**************************************************************************** * globale Variable ****************************************************************************/ /*-------------------------------------------------------------------------- * Files *-------------------------------------------------------------------------*/ FILE *fdyn,*fwlo,*f1y,*f3y,*fo4; /* ** fdyn= z88.dyn ** fwlo= z88e.log ** f1y= z88o1.bny ** f3y= z88o3.bny ** fo4= z88o4.txt */ char cdyn[8] = "z88.dyn"; char clog[9] = "z88e.log"; char c1y[10] = "z88o1.bny"; char c3y[10] = "z88o3.bny"; char co4[10] = "z88o4.txt"; /*-------------------------------------------------------------------------- * Pointer *-------------------------------------------------------------------------*/ FR_DOUBLEAY se; FR_DOUBLEAY u; FR_DOUBLEAY x; FR_DOUBLEAY y; FR_DOUBLEAY z; FR_DOUBLEAY emod; FR_DOUBLEAY rnue; FR_DOUBLEAY qpara; FR_DOUBLEAY riyy; FR_DOUBLEAY eyy; FR_DOUBLEAY rizz; FR_DOUBLEAY ezz; FR_DOUBLEAY rit; FR_DOUBLEAY wt; FR_DOUBLEAY fsum1; FR_DOUBLEAY fsum2; FR_DOUBLEAY fsum3; FR_DOUBLEAY fsum4; FR_DOUBLEAY fsum5; FR_DOUBLEAY fsum6; FR_INT4AY koi; FR_INT4AY ifrei; FR_INT4AY ioffs; FR_INT4AY koffs; FR_INT4AY ityp; FR_INT4AY ivon; FR_INT4AY ibis; FR_INT4AY intord; /*-------------------------------------------------------------------------- * Arrays *-------------------------------------------------------------------------*/ FR_DOUBLE f[61]; /* 61 ist MAXFE */ FR_DOUBLE ul[61]; /* 61 ist MAXFE */ FR_INT4 mcomp[21]; /* 21 ist MAXPA */ /* Diese Arrays werden in HEXA88,LQUA88,QSHE88 und CSHE88 verwendet */ FR_DOUBLE xk[21], yk[21], zk[21]; /* 21 ist MAXPA , HEXA88 */ FR_DOUBLE h[21]; /* 21 ist MAXPA , HEXA88 */ FR_DOUBLE b[361]; /* ist 6 x 60 +1, HEXA88 */ FR_DOUBLE xx[61]; /* ist 3 x 20 +1, HEXA88 */ FR_DOUBLE d[37]; /* ist 6 x 6 +1, HEXA88 */ FR_DOUBLE p[61]; /* ist 3 x 20 +1, HEXA88 */ /* fuer Plattenberechnung */ FR_DOUBLE be[49]; /* fuer 16-Knoten Platte */ FR_DOUBLE hi[49]; FR_DOUBLE bbi[145]; FR_DOUBLE bsv[97]; FR_DOUBLE dbi[10]; FR_DOUBLE dsv[5]; /*-------------------------------------------------------------------------- * Variable *-------------------------------------------------------------------------*/ FR_DOUBLE emode,rnuee,qparae,riyye,eyye,rizze,ezze,rite,wte; FR_INT4 intore,nel,ktyp; FR_INT4 IDYNMEM, LANG; FR_INT4 ndim,nkp,ne,nfg,neg,nfgp1,nkoi,ibflag,ipflag; FR_INT4 mxknot,mxfrei,mxfe; /*-------------------------------------------------------------------------- * vorbelegte Variable *-------------------------------------------------------------------------*/ FR_INT4 MAXNFG=0,MAXK=0,MAXE=0,MAXKOI=0,MAXNEG=0; FR_INT4 MAXESM=3600; /*********************************************************************** * hier beginnt z88e ***********************************************************************/ int main(void) { FR_INT4 i,k,j; int iret; /*--------------------------------------------------------------------------- * Die Sprache feststellen *--------------------------------------------------------------------------*/ LANG = 0; iret= lan88e(); if(iret != 0) { ale88e(iret); stop88e(); } /*--------------------------------------------------------------------------- * Ueberschrift *--------------------------------------------------------------------------*/ who88e(); /*--------------------------------------------------------------------------- * dyn88e starten: anfordern dynamic memory *--------------------------------------------------------------------------*/ iret= dyn88e(); if(iret != 0) { ale88e(iret); stop88e(); } /*--------------------------------------------------------------------------- * File Z88O4.TXT oeffnen *--------------------------------------------------------------------------*/ fo4= fopen(co4,"w+"); if(fo4 == NULL) { wlog88e(0,LOG_NOO4); fclose(fwlo); ale88e(AL_NOO4); stop88e(); } rewind(fo4); /*--------------------------------------------------------------------------- * riy88 liest Z88O1.BNY und Z88O3.BNY ein *--------------------------------------------------------------------------*/ wlog88e(0,LOG_RIY88); iret= riy88(); if(iret != 0) { ale88e(iret); stop88e(); } /*********************************************************************** * Berechnen der Elementsteifigkeitsmatrizen, Compilation ***********************************************************************/ wrim88e(0,TX_KNOTENK); wlog88e(0,LOG_KNOTENK); for(i= 1;i <= nkp;i++) /* schadet nicht */ { fsum1[i]= 0.; fsum2[i]= 0.; fsum3[i]= 0.; fsum4[i]= 0.; fsum5[i]= 0.; fsum6[i]= 0.; } if(LANG == 1) { fprintf(fo4,"Ausgabedatei Z88O4.TXT : Knotenkraefte, erzeugt mit Z88E V10\n"); fprintf(fo4," *************\n\n"); fprintf(fo4,"Zunaechst die Knotenkraefte elementweise\n"); fprintf(fo4,"----------------------------------------\n"); } if(LANG == 2) { fprintf(fo4,"output file Z88O4.TXT : nodal forces, computed by Z88E V10\n"); fprintf(fo4," ************\n\n"); fprintf(fo4,"starting with nodal forces computed for each element\n"); fprintf(fo4,"----------------------------------------------------\n"); } for(k= 1;k <= ne;k++) { /*---------------------------------------------------------------------- * es ist ityp(k): 1 isoparametrischer hexaeder, linearer ansatz * 2 balken * 3 scheibe, dreieckig, quadratischer ansatz * 4 raeumlicher stab * 5 welle * 6 torus, dreieckig, linearer ansatz * 7 isoparametrische scheibe, quadratischer ansatz * 8 isoparametrischer torus , quadratischer ansatz * 9 ebenes stabelement * 10 isop. serendipity hexaeder, quadratischer ans. * 11 isoparametrische scheibe, kubischer ansatz * 12 isoparametrischer torus, kubischer ansatz * 13 ebenes balkenelement * 14 krummlinige dreiecksscheibe, quadr. ansatz * 15 krummliniger dreieckstorus, quadr. ansatz * 16 krummliniger tetraeder, quadr. ansatz * 17 tetraeder, linearer ansatz * 18 6-knoten serendipity reissner mindlin platte * 19 16-knoten lagrange reissner mindlin platte * 20 8-knoten serendipity reissner mindlin platte *---------------------------------------------------------------------*/ /*---------------------------------------------------------------------- * zunaechst E-modul und Nue feststellen *---------------------------------------------------------------------*/ for(i = 1;i <= neg;i++) { if(k >= ivon[i] && k <= ibis[i]) { emode= emod[i]; rnuee= rnue[i]; intore=intord[i]; qparae= qpara[i]; if(ibflag == 1) { riyye= riyy[i]; eyye = eyy[i]; rizze= rizz[i]; ezze = ezz[i]; rite = rit[i]; wte = wt[i]; } if(ipflag != 0) { riyye= riyy[i]; } goto L60; } } L60:; /*---------------------------------------------------------------------- * Start 8-Knoten Quader, isoparametrischer linearer Ansatz *---------------------------------------------------------------------*/ if(ityp[k]== 1) { wtyp88e(k,1); /*---------------------------------------------------------------------- * isopara 8-Knoten Hexaeder: zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ for(i = 1;i <= 8;i++) { xk[i] = x [koi[koffs[k]+i-1]]; yk[i] = y [koi[koffs[k]+i-1]]; zk[i] = z [koi[koffs[k]+i-1]]; } /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer iso-8-Knoten Hexaeder berechen *---------------------------------------------------------------------*/ iret= lqua88(); if(iret != 0) { wlog88e(0,LOG_JACNEG); fclose(fwlo); ale88e(iret); stop88e(); } /*---------------------------------------------------------------------- * Compilation fuer lqua88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ for(i = 1;i <= 8;i++) mcomp[i]= ioffs[ koi[koffs[k]+i-1]] -1; mxknot= 8; mxfrei= 3; mxfe = 24; forc88(); if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 8-K Hexaeder",k); fprintf(fo4,"\nKnoten F(1) F(2) F(3)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 8-n hexahedron",k); fprintf(fo4,"\nnode F(1) F(2) F(3)"); } j= 1; for(i = 1;i <= 8;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1],f[j+2]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; fsum3[koi[koffs[k]+i-1]] += f[j+2]; j+= 3; } goto L7000; /*---------------------------------------------------------------------- * Ende isopara-8-Knoten Hexaeder *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start Balkenelement *---------------------------------------------------------------------*/ else if(ityp[k]== 2) { wtyp88e(k,2); /*---------------------------------------------------------------------- * Balkenelement : zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ xk[1] = x [koi[koffs[k]]]; yk[1] = y [koi[koffs[k]]]; zk[1] = z [koi[koffs[k]]]; xk[2] = x [koi[koffs[k]+1]]; yk[2] = y [koi[koffs[k]+1]]; zk[2] = z [koi[koffs[k]+1]]; /*---------------------------------------------------------------------- * Elementsteifigkeitsmatrix fuer Balkenelement berechen *---------------------------------------------------------------------*/ balk88(); /*---------------------------------------------------------------------- * Compilation fuer balk88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ mcomp[1]= ioffs[koi[koffs[k] ]] -1; mcomp[2]= ioffs[koi[koffs[k]+1]] -1; mxknot= 2; mxfrei= 6; mxfe = 12; forc88(); if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = Balken im Raum",k); fprintf(fo4,"\nKnoten F(1) F(2) F(3)\ F(4) F(5) F(6)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = beam in space",k); fprintf(fo4,"\nnode F(1) F(2) F(3)\ F(4) F(5) F(6)"); } j= 1; for(i = 1;i <= 2;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE %+#13.5lE\ %+#13.5lE %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1],f[j+2],f[j+3],f[j+4],f[j+5]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; fsum3[koi[koffs[k]+i-1]] += f[j+2]; fsum4[koi[koffs[k]+i-1]] += f[j+3]; fsum5[koi[koffs[k]+i-1]] += f[j+4]; fsum6[koi[koffs[k]+i-1]] += f[j+5]; j+= 6; } goto L7000; /*---------------------------------------------------------------------- * Ende Balkenelement *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start Dreieck-Scheibenelement, quadratischer Ansatz *---------------------------------------------------------------------*/ else if(ityp[k]== 3) { wtyp88e(k,3); /*---------------------------------------------------------------------- * Scheibenelement: zutreffende Koordinaten bestimmen * die Mittenknoten werden nicht gebraucht *---------------------------------------------------------------------*/ xk[1] = x [koi[koffs[k] ]]; yk[1] = y [koi[koffs[k] ]]; xk[2] = x [koi[koffs[k]+1]]; yk[2] = y [koi[koffs[k]+1]]; xk[3] = x [koi[koffs[k]+2]]; yk[3] = y [koi[koffs[k]+2]]; /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer Scheibenelement berechnen *---------------------------------------------------------------------*/ iret= shei88(); if(iret != 0) { wlog88e(0,LOG_JACNEG); fclose(fwlo); ale88e(iret); stop88e(); } /*----------------------------------------------------------------------- * Compilation fuer shei88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ for(i = 1;i <= 6;i++) mcomp[i]= ioffs[ koi[koffs[k]+i-1]] -1; mxknot= 6; mxfrei= 2; mxfe = 12; forc88(); if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 6-K Scheibe",k); fprintf(fo4,"\nKnoten F(1) F(2)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 6-n plane stress",k); fprintf(fo4,"\nnode F(1) F(2)"); } j= 1; for(i = 1;i <= 6;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; j+= 2; } goto L7000; /*---------------------------------------------------------------------- * Ende 6-k Scheibenelement *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start raeumliches Stabelement *---------------------------------------------------------------------*/ else if(ityp[k]== 4) { wtyp88e(k,4); /*---------------------------------------------------------------------- * raeumliches Stabelement : zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ xk[1] = x [koi[koffs[k]]]; yk[1] = y [koi[koffs[k]]]; zk[1] = z [koi[koffs[k]]]; xk[2] = x [koi[koffs[k]+1]]; yk[2] = y [koi[koffs[k]+1]]; zk[2] = z [koi[koffs[k]+1]]; /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer raeumliches Stabelement berechen *---------------------------------------------------------------------*/ stab88(); /*---------------------------------------------------------------------- * Compilation fuer stab88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ mcomp[1]= ioffs[koi[koffs[k] ]] -1; mcomp[2]= ioffs[koi[koffs[k]+1]] -1; mxknot= 2; mxfrei= 3; mxfe = 6; forc88(); if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = Stab im Raum",k); fprintf(fo4,"\nKnoten F(1) F(2) F(3)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = truss in space",k); fprintf(fo4,"\nnode F(1) F(2) F(3)"); } j= 1; for(i = 1;i <= 2;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1],f[j+2]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; fsum3[koi[koffs[k]+i-1]] += f[j+2]; j+= 3; } goto L7000; /*---------------------------------------------------------------------- * Ende raeumliches Stabelement *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start Wellenelement *---------------------------------------------------------------------*/ else if(ityp[k]== 5) { wtyp88e(k,5); /*---------------------------------------------------------------------- * Wellenelement : zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ xk[1] = x [koi[koffs[k]]]; xk[2] = x [koi[koffs[k]+1]]; /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer Wellenelement berechen *---------------------------------------------------------------------*/ well88(); /*---------------------------------------------------------------------- * Compilation fuer well88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ mcomp[1]= ioffs[ koi[koffs[k] ]] -1; mcomp[2]= ioffs[ koi[koffs[k]+1]] -1; mxknot= 2; mxfrei= 6; mxfe = 12; forc88(); if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = Welle",k); fprintf(fo4,"\nKnoten F(1) F(2) F(3)\ F(4) F(5) F(6)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = shaft",k); fprintf(fo4,"\nnode F(1) F(2) F(3)\ F(4) F(5) F(6)"); } j= 1; for(i = 1;i <= 2;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE %+#13.5lE\ %+#13.5lE %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1],f[j+2],f[j+3],f[j+4],f[j+5]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; fsum3[koi[koffs[k]+i-1]] += f[j+2]; fsum4[koi[koffs[k]+i-1]] += f[j+3]; fsum5[koi[koffs[k]+i-1]] += f[j+4]; fsum6[koi[koffs[k]+i-1]] += f[j+5]; j+= 6; } goto L7000; /*---------------------------------------------------------------------- * Ende Wellenelement *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start Toruselement, linearer Ansatz *---------------------------------------------------------------------*/ else if(ityp[k]== 6) { wtyp88e(k,6); /*---------------------------------------------------------------------- * Toruselement : zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ xk[1] = x [koi[koffs[k] ]]; yk[1] = y [koi[koffs[k] ]]; xk[2] = x [koi[koffs[k]+1]]; yk[2] = y [koi[koffs[k]+1]]; xk[3] = x [koi[koffs[k]+2]]; yk[3] = y [koi[koffs[k]+2]]; /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer Toruselement berechen *---------------------------------------------------------------------*/ ring88(); /*---------------------------------------------------------------------- * Compilation fuer ring88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ mcomp[1]= ioffs[koi[koffs[k] ]] -1; mcomp[2]= ioffs[koi[koffs[k]+1]] -1; mcomp[3]= ioffs[koi[koffs[k]+2]] -1; mxknot= 3; mxfrei= 2; mxfe = 6; forc88(); if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 3-K Torus",k); fprintf(fo4,"\nKnoten F(1) F(2)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 3-n torus",k); fprintf(fo4,"\nnode F(1) F(2)"); } j= 1; for(i = 1;i <= 3;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; j+= 2; } goto L7000; /*---------------------------------------------------------------------- * Ende Toruselement *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start 8-knoten Scheibe/Torus, isoparametrischer quadratischer Ansatz *---------------------------------------------------------------------*/ else if(ityp[k]== 7 || ityp[k]== 8) { if (ityp[k]== 7) { wtyp88e(k,7); ktyp= 2; } else { wtyp88e(k,8); ktyp= 0; } /*---------------------------------------------------------------------- * isopara 8-Knoten Scheibe/Torus: zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ for(i = 1;i <= 8;i++) { xk[i] = x [koi[koffs[k]+i-1]]; yk[i] = y [koi[koffs[k]+i-1]]; } /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer iso-8-knoten Scheibe/Tor berechen *---------------------------------------------------------------------*/ iret= qshe88(); if(iret != 0) { wlog88e(0,LOG_JACNEG); fclose(fwlo); ale88e(iret); stop88e(); } /*---------------------------------------------------------------------- * Compilation fuer qshe88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ for(i = 1;i <= 8;i++) mcomp[i]= ioffs[ koi[koffs[k]+i-1]] -1; mxknot= 8; mxfrei= 2; mxfe = 16; forc88(); if(ktyp == 2) { if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 8-K Scheibe",k); fprintf(fo4,"\nKnoten F(1) F(2)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 8-n plane stress",k); fprintf(fo4,"\nnode F(1) F(2)"); } } if(ktyp == 0) { if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 8-K Torus",k); fprintf(fo4,"\nKnoten F(1) F(2)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 8-n torus",k); fprintf(fo4,"\nnode F(1) F(2)"); } } j= 1; for(i = 1;i <= 8;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; j+= 2; } goto L7000; /*---------------------------------------------------------------------- * Ende isopara-8-Knoten Scheibe/Torus *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start ebenes Stabelement *---------------------------------------------------------------------*/ else if(ityp[k]== 9) { wtyp88e(k,9); /*---------------------------------------------------------------------- * ebenes Stabelement : zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ xk[1] = x [koi[koffs[k] ]]; yk[1] = y [koi[koffs[k] ]]; xk[2] = x [koi[koffs[k]+1]]; yk[2] = y [koi[koffs[k]+1]]; /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer e. Stabelement berechen *---------------------------------------------------------------------*/ esta88(); /*---------------------------------------------------------------------- * Compilation fuer esta88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ mcomp[1]= ioffs[ koi[koffs[k] ]] -1; mcomp[2]= ioffs[ koi[koffs[k]+1]] -1; mxknot= 2; mxfrei= 2; mxfe = 4; forc88(); if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = Stab in der Ebene",k); fprintf(fo4,"\nKnoten F(1) F(2)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = truss in plane",k); fprintf(fo4,"\nnode F(1) F(2)"); } j= 1; for(i = 1;i <= 2;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; j+= 2; } goto L7000; /*---------------------------------------------------------------------- * Ende ebenes Stabelement *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start 20-Knoten Hexaeder, isoparametrischer quadratischer Ansatz *---------------------------------------------------------------------*/ else if(ityp[k]== 10) { wtyp88e(k,10); /*---------------------------------------------------------------------- * isopara 20-Knoten Hexaeder: zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ for(i = 1;i <= 20;i++) { xk[i] = x [koi[koffs[k]+i-1]]; yk[i] = y [koi[koffs[k]+i-1]]; zk[i] = z [koi[koffs[k]+i-1]]; } /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer iso-20-Knoten Hexaeder berechen *---------------------------------------------------------------------*/ iret= hexa88(); if(iret != 0) { wlog88e(0,LOG_JACNEG); fclose(fwlo); ale88e(iret); stop88e(); } /*---------------------------------------------------------------------- * Compilation fuer hexa88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ for(i = 1;i <= 20;i++) mcomp[i]= ioffs[ koi[koffs[k]+i-1]] -1; mxknot= 20; mxfrei= 3; mxfe = 60; forc88(); if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 20-K Hexaeder",k); fprintf(fo4,"\nKnoten F(1) F(2) F(3)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 20-n hexahedron",k); fprintf(fo4,"\nnode F(1) F(2) F(3)"); } j= 1; for(i = 1;i <= 20;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1],f[j+2]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; fsum3[koi[koffs[k]+i-1]] += f[j+2]; j+= 3; } goto L7000; /*---------------------------------------------------------------------- * Ende isopara-20-Knoten Hexaeder *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start 12-knoten Scheibe/Torus, isoparametrischer kubischer Ansatz *---------------------------------------------------------------------*/ else if(ityp[k]== 11 || ityp[k]== 12) { if (ityp[k]== 11) { wtyp88e(k,11); ktyp= 2; } else { wtyp88e(k,12); ktyp= 0; } /*---------------------------------------------------------------------- * isopara 12-Knoten Scheibe/Torus: zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ for(i = 1;i <= 12;i++) { xk[i] = x [koi[koffs[k]+i-1]]; yk[i] = y [koi[koffs[k]+i-1]]; } /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer iso-12-Knoten Scheibe/Tor ber. *---------------------------------------------------------------------*/ iret= cshe88(); if(iret != 0) { wlog88e(0,LOG_JACNEG); fclose(fwlo); ale88e(iret); stop88e(); } /*---------------------------------------------------------------------- * Compilation fuer cshe88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ for(i = 1;i <= 12;i++) mcomp[i]= ioffs[ koi[koffs[k]+i-1]] -1; mxknot= 12; mxfrei= 2; mxfe = 24; forc88(); if(ktyp == 2) { if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 12-K Scheibe",k); fprintf(fo4,"\nKnoten F(1) F(2)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 12-n plane stress",k); fprintf(fo4,"\nnode F(1) F(2)"); } } if(ktyp == 0) { if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 12-K Torus",k); fprintf(fo4,"\nKnoten F(1) F(2)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 12-n torus",k); fprintf(fo4,"\nnode F(1) F(2)"); } } j= 1; for(i = 1;i <= 12;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; j+= 2; } goto L7000; /*---------------------------------------------------------------------- * Ende isopara-12-Knoten Scheibe/Torus *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start ebenes Balkenelement *---------------------------------------------------------------------*/ else if(ityp[k]== 13) { wtyp88e(k,13); /*---------------------------------------------------------------------- * Balkenelement : zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ xk[1] = x [koi[koffs[k]]]; yk[1] = y [koi[koffs[k]]]; xk[2] = x [koi[koffs[k]+1]]; yk[2] = y [koi[koffs[k]+1]]; /*---------------------------------------------------------------------- * Elementsteifigkeitsmatrix fuer ebenes Balkenelement berechen *---------------------------------------------------------------------*/ ebal88(); /*---------------------------------------------------------------------- * Compilation fuer ebal88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ mcomp[1]= ioffs[koi[koffs[k] ]] -1; mcomp[2]= ioffs[koi[koffs[k]+1]] -1; mxknot= 2; mxfrei= 3; mxfe = 6; forc88(); if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = Balken in der Ebene",k); fprintf(fo4,"\nKnoten F(1) F(2) F(3)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = beam in plane",k); fprintf(fo4,"\nnode F(1) F(2) F(3)"); } j= 1; for(i = 1;i <= 2;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1],f[j+2]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; fsum3[koi[koffs[k]+i-1]] += f[j+2]; j+= 3; } goto L7000; /*---------------------------------------------------------------------- * Ende ebenes Balkenelement *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start 6-knoten Dreieck-Scheibe/Torus, isop. quadratischer Ansatz *---------------------------------------------------------------------*/ else if(ityp[k]== 14 || ityp[k]== 15) { if (ityp[k]== 14) { wtyp88e(k,14); ktyp= 2; } else { wtyp88e(k,15); ktyp= 0; } /*---------------------------------------------------------------------- * iso 6-Knoten Dreieck-Scheibe/Torus: zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ for(i = 1;i <= 6;i++) { xk[i] = x [koi[koffs[k]+i-1]]; yk[i] = y [koi[koffs[k]+i-1]]; } /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer iso-6-knoten Scheibe/Tor berechnen *---------------------------------------------------------------------*/ iret= isod88(); if(iret != 0) { wlog88e(0,LOG_JACNEG); fclose(fwlo); ale88e(iret); stop88e(); } /*---------------------------------------------------------------------- * Compilation fuer isod88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ for(i = 1;i <= 6;i++) mcomp[i]= ioffs[ koi[koffs[k]+i-1]] -1; mxknot= 6; mxfrei= 2; mxfe = 12; forc88(); if(ktyp == 2) { if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 6-K isop. Scheibe",k); fprintf(fo4,"\nKnoten F(1) F(2)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 6-n isop. plane stress",k); fprintf(fo4,"\nnode F(1) F(2)"); } } if(ktyp == 0) { if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 6-K isop. Torus",k); fprintf(fo4,"\nKnoten F(1) F(2)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 6-n isop. torus",k); fprintf(fo4,"\nnode F(1) F(2)"); } } j= 1; for(i = 1;i <= 6;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; j+= 2; } goto L7000; /*---------------------------------------------------------------------- * Ende isopara-6-Knoten Dreieck-Scheibe/Torus *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start 10-Knoten Tetraeder, isoparametrischer quadratischer Ansatz *---------------------------------------------------------------------*/ else if(ityp[k]== 16) { wtyp88e(k,16); /*---------------------------------------------------------------------- * isopara 10-Knoten Tetraeder: zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ for(i = 1;i <= 10;i++) { xk[i] = x [koi[koffs[k]+i-1]]; yk[i] = y [koi[koffs[k]+i-1]]; zk[i] = z [koi[koffs[k]+i-1]]; } /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer iso-10-Knoten Tetraeder berechen *---------------------------------------------------------------------*/ iret= tetr88(); if(iret != 0) { wlog88e(0,LOG_JACNEG); fclose(fwlo); ale88e(iret); stop88e(); } /*---------------------------------------------------------------------- * Compilation fuer tetr88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ for(i = 1;i <= 10;i++) mcomp[i]= ioffs[ koi[koffs[k]+i-1]] -1; mxknot= 10; mxfrei= 3; mxfe = 30; forc88(); if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 10-K Tetraeder",k); fprintf(fo4,"\nKnoten F(1) F(2) F(3)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 10-n tetrahedron",k); fprintf(fo4,"\nnode F(1) F(2) F(3)"); } j= 1; for(i = 1;i <= 10;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1],f[j+2]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; fsum3[koi[koffs[k]+i-1]] += f[j+2]; j+= 3; } goto L7000; /*---------------------------------------------------------------------- * Ende isopara-10-Knoten Tetraeder *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start 4-Knoten Tetraeder, lineaer Ansatz *---------------------------------------------------------------------*/ else if(ityp[k]== 17) { wtyp88e(k,17); /*---------------------------------------------------------------------- * 4-Knoten Tetraeder: zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ for(i = 1;i <= 4;i++) { xk[i] = x [koi[koffs[k]+i-1]]; yk[i] = y [koi[koffs[k]+i-1]]; zk[i] = z [koi[koffs[k]+i-1]]; } /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer 4-Knoten Tetraeder berechen *---------------------------------------------------------------------*/ iret= spur88(); if(iret != 0) { wlog88e(0,LOG_JACNEG); fclose(fwlo); ale88e(iret); stop88e(); } /*---------------------------------------------------------------------- * Compilation fuer spur88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ for(i = 1;i <= 4;i++) mcomp[i]= ioffs[ koi[koffs[k]+i-1]] -1; mxknot= 4; mxfrei= 3; mxfe = 12; forc88(); if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 4-K Tetraeder",k); fprintf(fo4,"\nKnoten F(1) F(2) F(3)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 4-n tetrahedron",k); fprintf(fo4,"\nnode F(1) F(2) F(3)"); } j= 1; for(i = 1;i <= 4;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1],f[j+2]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; fsum3[koi[koffs[k]+i-1]] += f[j+2]; j+= 3; } goto L7000; /*---------------------------------------------------------------------- * Ende 4-Knoten Tetraeder *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start 6-Knoten Reissner- Mindlin Platte *---------------------------------------------------------------------*/ else if(ityp[k]== 18) { wtyp88e(k,18); /*---------------------------------------------------------------------- * 6-Knoten Plattenelement : zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ for(i = 1;i <= 6;i++) { xk[i] = x [koi[koffs[k]+i-1]]; yk[i] = y [koi[koffs[k]+i-1]]; } /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer 6-Knoten Plattenelement berechen *---------------------------------------------------------------------*/ iret= spla88(); if(iret != 0) { wlog88e(0,LOG_JACNEG); fclose(fwlo); ale88e(iret); stop88e(); } /*---------------------------------------------------------------------- * Compilation fuer spla88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ for(i = 1;i <= 6;i++) mcomp[i]= ioffs[ koi[koffs[k]+i-1]] -1; mxknot= 6; mxfrei= 3; mxfe = 18; forc88(); if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 6-K Mindlin- Platte",k); fprintf(fo4,"\nKnoten F(1) F(2) F(3)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 6-n Mindlin plate",k); fprintf(fo4,"\nnode F(1) F(2) F(3)"); } j= 1; for(i = 1;i <= 6;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1],f[j+2]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; fsum3[koi[koffs[k]+i-1]] += f[j+2]; j+= 3; } goto L7000; /*---------------------------------------------------------------------- * Ende 6-Knoten Reissner- Mindlin Platte *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start 16-Knoten Reissner Mindlin Platte, kubischer Lagrange Ansatz *---------------------------------------------------------------------*/ else if(ityp[k]== 19) { wtyp88e(k,19); /*---------------------------------------------------------------------- * 16-Knoten Platte: zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ for(i = 1;i <= 16;i++) { xk[i] = x [koi[koffs[k]+i-1]]; yk[i] = y [koi[koffs[k]+i-1]]; } /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer 16-Knoten Platte berechen *---------------------------------------------------------------------*/ iret= hpla88(); if(iret != 0) { wlog88e(0,LOG_JACNEG); fclose(fwlo); ale88e(iret); stop88e(); } /*---------------------------------------------------------------------- * Compilation fuer hpla88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ for(i = 1;i <= 16;i++) mcomp[i]= ioffs[ koi[koffs[k]+i-1]] -1; mxknot= 16; mxfrei= 3; mxfe = 48; forc88(); if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 16-K Mindlin Platte",k); fprintf(fo4,"\nKnoten F(1) F(2) F(3)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 16-n Mindlin plate",k); fprintf(fo4,"\nnode F(1) F(2) F(3)"); } j= 1; for(i = 1;i <= 16;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1],f[j+2]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; fsum3[koi[koffs[k]+i-1]] += f[j+2]; j+= 3; } goto L7000; /*---------------------------------------------------------------------- * Ende 16-Knoten Lagrange Reissner Mindlin Platte *---------------------------------------------------------------------*/ } /*---------------------------------------------------------------------- * Start 8-Knoten Reissner- Mindlin Serendipity Platte *---------------------------------------------------------------------*/ else if(ityp[k]== 20) { wtyp88e(k,20); /*---------------------------------------------------------------------- * 8-Knoten Plattenelement : zutreffende Koordinaten bestimmen *---------------------------------------------------------------------*/ for(i = 1;i <= 8;i++) { xk[i] = x [koi[koffs[k]+i-1]]; yk[i] = y [koi[koffs[k]+i-1]]; } /*---------------------------------------------------------------------- * nun Elementsteifigkeitsmatrix fuer 8-Knoten Plattenelement berechen *---------------------------------------------------------------------*/ iret= apla88(); if(iret != 0) { wlog88e(0,LOG_JACNEG); fclose(fwlo); ale88e(iret); stop88e(); } /*---------------------------------------------------------------------- * Compilation fuer apla88, kompakte Speicherung mit Pointervektor *---------------------------------------------------------------------*/ for(i = 1;i <= 8;i++) mcomp[i]= ioffs[ koi[koffs[k]+i-1]] -1; mxknot= 8; mxfrei= 3; mxfe = 24; forc88(); if(LANG == 1) { fprintf(fo4,"\nElement # = %5ld Typ = 8-K Mindlin- Platte",k); fprintf(fo4,"\nKnoten F(1) F(2) F(3)"); } if(LANG == 2) { fprintf(fo4,"\nelement # = %5ld type = 8-n Mindlin plate",k); fprintf(fo4,"\nnode F(1) F(2) F(3)"); } j= 1; for(i = 1;i <= 8;i++) { fprintf(fo4,"\n%5ld %+#13.5lE %+#13.5lE %+#13.5lE", koi[koffs[k]+i-1],f[j],f[j+1],f[j+2]); fsum1[koi[koffs[k]+i-1]] += f[j]; fsum2[koi[koffs[k]+i-1]] += f[j+1]; fsum3[koi[koffs[k]+i-1]] += f[j+2]; j+= 3; } goto L7000; /*---------------------------------------------------------------------- * Ende 8-Knoten Reissner- Mindlin Serendipity Platte *---------------------------------------------------------------------*/ } L7000:; fprintf(fo4,"\n"); /*********************************************************************** * Ende der Schleife ueber alle Elemente ***********************************************************************/ } /* e 8000 */ /*********************************************************************** * Schleife ueber alle Knoten ***********************************************************************/ fprintf(fo4,"\n"); if(LANG == 1) { fprintf(fo4,"Nun die aufsummierten Knotenkraefte je Knoten\n"); fprintf(fo4,"---------------------------------------------\n\n"); fprintf(fo4,"Knoten F(1) F(2) F(3)\ F(4) F(5) F(6)\n\n"); } if(LANG == 2) { fprintf(fo4,"now the nodal sums for each node\n"); fprintf(fo4,"--------------------------------\n"); fprintf(fo4," node F(1) F(2) F(3)\ F(4) F(5) F(6)\n\n"); } for(i = 1;i <= nkp;i++) { fprintf(fo4, "%5ld %+#13.5lE %+#13.5lE %+#13.5lE %+#13.5lE %+#13.5lE %+#13.5lE\n", i,fsum1[i],fsum2[i],fsum3[i],fsum4[i],fsum5[i],fsum6[i]); } /*---------------------------------------------------------------------- * Ende Z88E *---------------------------------------------------------------------*/ fclose(fo4); wrim88e(0,TX_EXITZ88E); wlog88e(0,LOG_EXITZ88E); return(0); }