1 : C BBST2 SOURCE CHAT 05/01/12 21:36:37 5004 2 : SUBROUTINE BBST2(TRACE,NBPGAU,IFOUR,MELE,POIGAU,QSIGAU, 3 : & ETAGAU,DZEGAU,SHPTOT,NBNO,SHP,XE) 4 : c 5 : c 6 : c projection de la trace de epsi2 sur la base reduite des elements icq 7 : c 8 : c entree 9 : c trace : tableau contenant la valeur de trace de epsi2 en 10 : c chaque point de gauss 11 : c nbpgau : nombre de point de gauss 12 : c ifour : variable du ccoptio 13 : c mele : numero de l'element 14 : c poigau : poids d'integration 15 : c qsigau : 1ere composante sur l'element de reference 16 : c etagau : 2eme composante sur l'element de reference 17 : c dzegau : 3eme composante sur l'element de reference 18 : c shptot : valeurs des fonctions de forme 19 : c 20 : c sortie 21 : c trac : contient la correction à apporter à chaque composante diagonale de 22 : c la deformation 23 : c 24 : c 25 : IMPLICIT INTEGER(I-N) 26 : IMPLICIT REAL*8(A-H,O-Z) 27 : c 28 : DIMENSION TRACE(*),POIGAU(*),QSIGAU(*),ETAGAU(*),DZEGAU(*) 29 : & ,SHPTOT(6,NBNO,*),SHP(6,*),XE(3,*) 30 : DIMENSION P(3,3),T(3) 31 : 32 : c 33 : c elements icq4 et ict3 et ict6 34 : c 35 : IF (MELE.EQ.69.OR.MELE.EQ.70.OR.MELE.EQ.71) THEN 36 : TRINTG = 0.D0 37 : aire = 0.d0 38 : DO 10 I=1,NBPGAU 39 : DO 5 J=1,NBNO 40 : DO 5 K =1,6 41 : SHP(K,J)=SHPTOT(K,J,I) 42 : 5 CONTINUE 42 : CALL JACOBI(XE,SHP,2,NBNO,DJAC) 43 : TRINTG = TRINTG + POIGAU(I)*TRACE(I)*DJAC 44 : AIRE = AIRE +POIGAU(I)*DJAC 45 : 10 CONTINUE 46 : trintg = TRINTG / AIRE 47 : c 48 : COEF = 2.D0 49 : IF (IFOUR .EQ. 0 ) COEF = 3.D0 50 : DO 20 I=1,NBPGAU 51 : TRACE(I) = (TRINTG - TRACE(I))/COEF 52 : 20 CONTINUE 53 : RETURN 54 : c 55 : c élement ICQ8 56 : c 57 : ELSE IF ( MELE.EQ.72) THEN 58 : c 59 : c matrice P calculée analytiquement voir routine BBCALC 60 : c 61 : P(1,1) = 2.D0/3.D0*(XE(1,1)*(XE(2,2)-XE(2,8)) 62 : 1 +XE(1,2)*(XE(2,3)-XE(2,1)) 63 : 1 +XE(1,3)*(XE(2,4)-XE(2,2))+XE(1,4)*(XE(2,5)-XE(2,3)) 64 : 2 +XE(1,5)*(XE(2,6)-XE(2,4))+XE(1,6)*(XE(2,7)-XE(2,5)) 65 : 3 +XE(1,7)*(XE(2,8)-XE(2,6))+XE(1,8)*(XE(2,1)-XE(2,7))) 66 : 4 +1.D0/6.D0*((XE(1,1)-XE(1,5))*(XE(2,7)-XE(2,3)) 67 : 5 +(XE(1,3)-XE(1,7))*(XE(2,1)-XE(2,5))) 68 : 69 : P(1,2) = 5.D0/9.D0*(-XE(2,2)*(XE(1,1)+XE(1,3)) 70 : 1 +XE(1,2)*(XE(2,1)+XE(2,3))+XE(1,6)*(-XE(2,5)-XE(2,7)) 71 : 1 +XE(2,6)*(XE(1,5)+XE(1,7))) 72 : 2 +4.D0/9.D0*(XE(1,8)*(XE(2,7)+XE(2,2)-XE(2,6)-XE(2,1)) 73 : 2 +XE(2,8)*(-XE(1,7)-XE(1,2)+XE(1,6)+XE(1,1)) 74 : 2 +XE(1,4)*(XE(2,2)-XE(2,6)+XE(2,5)-XE(2,3)) 75 : 2 +XE(2,4)*(-XE(1,2)+XE(1,6)-XE(1,5)+XE(1,3))) 76 : 3 +7.D0/45.D0*(XE(1,2)*(XE(2,5)+XE(2,7)) 77 : 3 -XE(2,2)*(XE(1,5)+XE(1,7))-XE(1,6)*(XE(2,1)+XE(2,3)) 78 : 3 +XE(2,6)*(XE(1,1)+XE(1,3))) 79 : 4 +8.D0/15.D0*(XE(1,6)*XE(2,2)-XE(1,2)*XE(2,6)) 80 : 5 +7.D0/90.D0*(XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5) 81 : 5 +XE(1,5)*XE(2,3)-XE(1,1)*XE(2,7)) 82 : 6 +1.D0/30.D0*(-XE(1,7)*XE(2,3)+XE(1,3)*XE(2,7) 83 : 6 -XE(1,5)*XE(2,1)+XE(1,1)*XE(2,5)) 84 : 85 : 86 : P(1,3) = 5.D0/9.D0*(-XE(1,8)*(XE(2,7)+XE(2,1)) 87 : 1 +XE(2,8)*(XE(1,7)+XE(1,1)) 88 : 1 +XE(1,4)*(XE(2,3)+XE(2,5)) 89 : 1 -XE(2,4)*(XE(1,3)+XE(1,5))) 90 : 2 +4.D0/9.D0*((XE(1,8)-XE(1,4))*(XE(2,2)+XE(2,6)) 91 : 2 +(-XE(2,8)+XE(2,4))*(XE(1,2)+XE(1,6)) 92 : 2 +XE(1,2)*(XE(2,1)-XE(2,3))+XE(2,2)*(-XE(1,1)+XE(1,3)) 93 : 2 +XE(1,6)*(XE(2,7)-XE(2,5))+XE(2,6)*(-XE(1,7)+XE(1,5))) 94 : 3 +7.D0/45.D0*(-XE(1,8)*(XE(2,3)+XE(2,5)) 95 : 3 +XE(2,8)*(XE(1,3)+XE(1,5)) 96 : 3 +XE(1,4)*(XE(2,1)+XE(2,7)) 97 : 3 -XE(2,4)*(XE(1,1)+XE(1,7))) 98 : 4 +1.D0/30.D0*(-XE(1,7)*XE(2,3)+XE(1,3)*XE(2,7) 99 : 4 +XE(1,5)*XE(2,1)-XE(1,1)*XE(2,5)) 100 : 5 +7.D0/90.D0*(XE(1,7)*XE(2,5)-XE(1,3)*XE(2,1) 101 : 5 -XE(1,5)*XE(2,7)+XE(1,1)*XE(2,3)) 102 : 6 +8.D0/15.D0*(XE(1,8)*XE(2,4)-XE(1,4)*XE(2,8)) 103 : 104 : P(2,1) = P(1,2) 105 : 106 : P(2,2) = 16.D0/45.D0*(XE(2,6)*(XE(1,5)-XE(1,7)) 107 : 1 +XE(1,6)*(XE(2,7)-XE(2,5))+XE(1,2)*(XE(2,3)-XE(2,1)) 108 : 1 +XE(2,2)*(XE(1,1)-XE(1,3))) 109 : 2 +14.D0/45.D0*(XE(2,4)*(XE(1,3)-XE(1,5)) 110 : 2 +XE(1,4)*(XE(2,5)-XE(2,3))+XE(1,8)*(XE(2,1)-XE(2,7)) 111 : 2 +XE(2,8)*(XE(1,7)-XE(1,1))) 112 : 3 +8.D0/45.D0*(XE(1,4)*(XE(2,2)-XE(2,6)) 113 : 3 +XE(1,8)*(XE(2,6)-XE(2,2))+XE(1,2)*(XE(2,8)-XE(2,4)) 114 : 3 +XE(1,6)*(XE(2,4)-XE(2,8))) 115 : 4 +2.D0/45.D0*(XE(2,2)*(XE(1,5)-XE(1,7)) 116 : 4 +XE(1,2)*(XE(2,7)-XE(2,5)) +XE(1,6)*(XE(2,3)-XE(2,1)) 117 : 4 +XE(2,6)*(XE(1,1)-XE(1,3))) 118 : 5 +4.D0/45.D0*(XE(2,4)*(XE(1,1)-XE(1,7)) 119 : 5 +XE(1,4)*(XE(2,7)-XE(2,1)) +XE(1,8)*(XE(2,3)-XE(2,5)) 120 : 5 +XE(2,8)*(XE(1,5)-XE(1,3))) 121 : 6 +17.D0/90.D0*(XE(1,7)*XE(2,5)+XE(1,3)*XE(2,1) 122 : 6 -XE(1,5)*XE(2,7)-XE(1,1)*XE(2,3)) 123 : 7 +1.D0/90.D0*(-XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5) 124 : 7 +XE(1,5)*XE(2,3)+XE(1,1)*XE(2,7)) 125 : 126 : 127 : P(2,3) = 1.D0/3.D0*(XE(1,5)*(XE(2,6)-XE(2,4)) 128 : 1 +XE(1,8)*(XE(2,7)+XE(2,1))+XE(1,7)*(XE(2,6)-XE(2,8)) 129 : 1 -XE(1,2)*(XE(2,1)+XE(2,3))+XE(1,1)*(XE(2,2)-XE(2,8)) 130 : 1 +XE(1,4)*(XE(2,5)+XE(2,3))-XE(1,6)*(XE(2,7)+XE(2,5)) 131 : 1 +XE(1,3)*(XE(2,2)-XE(2,4))) 132 : 2 +4.D0/9.D0*(-XE(1,4)*(XE(2,2)+XE(2,6)) 133 : 2 -XE(1,8)*(XE(2,6)+XE(2,2))+XE(1,2)*(XE(2,8)+XE(2,4)) 134 : 2 +XE(1,6)*(XE(2,4)+XE(2,8))) 135 : 3 +1.D0/9.D0*(XE(1,7)*(XE(2,2)-XE(2,4)) 136 : 3 +XE(1,8)*(XE(2,3)+XE(2,5)) 137 : 3 +XE(1,5)*(XE(2,2)-XE(2,8))+XE(1,3)*(XE(2,6)-XE(2,8)) 138 : 3 -XE(1,2)*(XE(2,5)+XE(2,7))+XE(1,1)*(XE(2,6)-XE(2,4)) 139 : 3 +XE(1,4)*(XE(2,1)+XE(2,7))-XE(1,6)*(XE(2,1)+XE(2,3))) 140 : 141 : 142 : P(3,1) = P(1,3) 143 : 144 : P(3,2) = P(2,3) 145 : 146 : P(3,3) = 14.D0/45.D0*(XE(1,6)*(XE(2,7)-XE(2,5)) 147 : 1 +XE(2,6)*(XE(1,5)-XE(1,7))+XE(1,2)*(XE(2,3)-XE(2,1)) 148 : 1 +XE(2,2)*(XE(1,1)-XE(1,3))) 149 : 2 +16.D0/45.D0*(XE(1,8)*(XE(2,1)-XE(2,7)) 150 : 2 +XE(2,8)*(XE(1,7)-XE(1,1))+XE(1,4)*(XE(2,5)-XE(2,3)) 151 : 2 +XE(2,4)*(XE(1,3)-XE(1,5))) 152 : 3 +8.D0/45.D0*(XE(1,4)*(XE(2,2)-XE(2,6)) 153 : 3 +XE(1,8)*(XE(2,6)-XE(2,2))+XE(1,2)*(XE(2,8)-XE(2,4)) 154 : 3 +XE(1,6)*(XE(2,4)-XE(2,8))) 155 : 4 +2.D0/45.D0*(XE(2,4)*(XE(1,7)-XE(1,1)) 156 : 4 +XE(1,8)*(XE(2,5)-XE(2,3))+XE(2,8)*(-XE(1,5)+XE(1,3)) 157 : 4 +XE(1,4)*(XE(2,1)-XE(2,7))) 158 : 5 +4.D0/45.D0*(XE(2,2)*(XE(1,7)-XE(1,5)) 159 : 5 +XE(1,2)*(XE(2,5)-XE(2,7))+XE(1,6)*(XE(2,1)-XE(2,3)) 160 : 5 +XE(2,6)*(XE(1,3)-XE(1,1))) 161 : 6 +1.D0/90.D0*(-XE(1,5)*XE(2,7)+XE(1,3)*XE(2,1) 162 : 6 +XE(1,7)*XE(2,5)-XE(1,1)*XE(2,3)) 163 : 7 +17.D0/90.D0*(-XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5) 164 : 7 +XE(1,5)*XE(2,3)+XE(1,1)*XE(2,7)) 165 : 166 : c 167 : c calcul des produits, , 168 : c 169 : T(1) = 0.D0 170 : T(2) = 0.D0 171 : T(3) = 0.D0 172 : DO 40 I=1,NBPGAU 173 : DO 30 J=1,NBNO 174 : DO 30 K =1,6 175 : SHP(K,J)=SHPTOT(K,J,I) 176 : 30 CONTINUE 176 : CALL JACOBI(XE,SHP,2,NBNO,DJAC) 177 : T(1) = T(1) + POIGAU(I)*TRACE(I)*DJAC 178 : T(2) = T(2) + POIGAU(I)*TRACE(I)*DJAC*QSIGAU(I) 179 : T(3) = T(3) + POIGAU(I)*TRACE(I)*DJAC*ETAGAU(I) 180 : 40 CONTINUE 181 : * 182 : c resolution du systme P.X=T 183 : c 183 : CALL GAUSSJ(P,3,3,T,1,1) 184 : * 185 : c caclul du terme correctif 186 : c 187 : COEF = 2.D0 188 : IF (IFOUR .EQ. 0 ) COEF = 3.D0 189 : DO 50 I=1,NBPGAU 190 : TRACE(I)=(T(1)+QSIGAU(I)*T(2)+ETAGAU(I)*T(3)-TRACE(I))/COEF 191 : 50 CONTINUE 192 : c 193 : RETURN 194 : c 195 : ENDIF 196 : c 197 : c pas de correction pour les autres éléments 198 : c 199 : DO 70 I=1,NBPGAU 200 : TRACE(I)=0.D0 201 : 70 CONTINUE 202 : c 203 : c 204 : RETURN 205 : END 206 : 207 : © Cast3M 2003 - All rights reserved.
Disclaimer