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