1 : C CUBTWE SOURCE CHAT 05/01/12 22:32:47 5004 2 : SUBROUTINE CUBTWE(ICM,XE,XEC,XX,SHPC,SGA,QSI,IG) 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8(A-H,O-Z) 5 : DIMENSION A(3,20),A1(3,8),A2(2,8),Q(3,20),P(3,8), 6 : & AJ(3,3),X(3),XE(3,*),XEC(3,2),XX(3,2),AL(3),DS(6),XT(3), 7 : & BJ(3,3),F(3),SHPC(6,20,*),SGA(3,5,5),QSI(3,2),DT(3),SIB(3) 8 : DATA EPSI/1.D-7/, UN/1.D0/,TWO/2.D0/ 9 : DATA DS/-1.D0, 1.D0, -1.D0, 1.D0, -1.D0,-1.D0/ 10 : C --- ORDRE: 1 3 5 7 13 15 17 19; 2 6 14 18; 4 8 16 20; 9 10 11 12 11 : DATA Q/-1.D0,-1.D0,-1.D0, 1.D0,-1.D0,-1.D0, 1.D0, 1.D0,-1.D0, 12 : & -1.D0, 1.D0,-1.D0, -1.D0,-1.D0, 1.D0, 1.D0,-1.D0, 1.D0, 13 : & 1.D0, 1.D0, 1.D0, -1.D0, 1.D0, 1.D0, 0.D0,-1.D0,-1.D0, 14 : & 0.D0, 1.D0,-1.D0, 0.D0,-1.D0, 1.D0, 0.D0, 1.D0, 1.D0, 15 : & 1.D0, 0.D0,-1.D0, -1.D0, 0.D0,-1.D0, 1.D0, 0.D0, 1.D0, 16 : & -1.D0, 0.D0, 1.D0, -1.D0,-1.D0, 0.D0, 1.D0,-1.D0, 0.D0, 17 : & 1.D0, 1.D0, 0.D0, -1.D0, 1.D0, 0.D0/ 18 : D12=0.D0 19 : DO 5 I=1,3 20 : D12=D12+(XEC(I,2)-XEC(I,1))**2 21 : DO 5 J=1,20 22 : 5 A(I,J)=0.D0 23 : D12=DSQRT(D12) 24 : DO 6 I=1,4 25 : K=2*I-1 26 : DO 6 J=1,3 27 : P(J,I)=XE(J,K)*0.125D0 28 : 6 P(J,I+4)=XE(J,K+12)*0.125D0 29 : C IG: NUMERO DU POINT DE GAUSS 30 : DO 10 J=1,8 31 : DO 10 I=1,3 32 : CO=P(I,J) 33 : DO 14 JH=1,3 34 : 14 X(JH)=CO*Q(JH,J) 35 : A(I,1)=A(I,1)-2.D0*CO 36 : A(I,2)=A(I,2)-X(1) 37 : A(I,3)=A(I,3)-X(2) 38 : A(I,4)=A(I,4)-X(3) 39 : A(I,5)=A(I,5)+CO 40 : A(I,7)=A(I,7)+CO 41 : A(I,9)=A(I,9)+CO 42 : A(I,11)=A(I,11)+X(2) 43 : A(I,12)=A(I,12)+X(1) 44 : A(I,13)=A(I,13)+X(3) 45 : A(I,14)=A(I,14)+X(2) 46 : A(I,15)=A(I,15)+X(1) 47 : A(I,16)=A(I,16)+X(3) 48 : A(I,17)=A(I,17)+X(1)*Q(2,J)*Q(3,J) 49 : A(I,18)=A(I,18)+X(2)*Q(3,J) 50 : A(I,19)=A(I,19)+X(1)*Q(3,J) 51 : 10 A(I,20)=A(I,20)+X(1)*Q(2,J) 52 : DO 11 I=1,3 53 : P(I,1)=XE(I,2)*0.25D0 54 : P(I,2)=XE(I,6)*0.25D0 55 : P(I,3)=XE(I,14)*0.25D0 56 : 11 P(I,4)=XE(I,18)*0.25D0 57 : DO 12 J=1,4 58 : K=J+8 59 : DO 12 I=1,3 60 : CO=P(I,J) 61 : CC=CO*Q(2,K) 62 : A(I,1)=A(I,1)+CO 63 : A(I,3)=A(I,3)+CC 64 : A(I,4)=A(I,4)+CO*Q(3,K) 65 : A(I,5)=A(I,5)-CO 66 : A(I, 8)=A(I, 8)+CC*Q(3,K) 67 : A(I,11)=A(I,11)-CC 68 : A(I,16)=A(I,16)-CO*Q(3,K) 69 : 12 A(I,18)=A(I,18)-CC*Q(3,K) 70 : DO 21 I=1,3 71 : P(I,1)=XE(I,4)*0.25D0 72 : P(I,2)=XE(I,8)*0.25D0 73 : P(I,3)=XE(I,16)*0.25D0 74 : 21 P(I,4)=XE(I,20)*0.25D0 75 : DO 22 J=1,4 76 : K=J+12 77 : DO 22 I=1,3 78 : CO=P(I,J) 79 : CC=CO*Q(1,K) 80 : A(I,1)=A(I,1)+CO 81 : A(I,2)=A(I,2)+CC 82 : A(I,4)=A(I,4)+CO*Q(3,K) 83 : A(I,7)=A(I,7)-CO 84 : A(I,10)=A(I,10)+CC*Q(3,K) 85 : A(I,12)=A(I,12)-CC 86 : A(I,13)=A(I,13)-CO*Q(3,K) 87 : 22 A(I,19)=A(I,19)-CC*Q(3,K) 88 : DO 31 I=1,3 89 : DO 31 J=1,4 90 : 31 P(I,J)=XE(I,J+8)*0.25D0 91 : DO 32 J=1,4 92 : K=J+16 93 : DO 32 I=1,3 94 : CO=P(I,J) 95 : CC=CO*Q(1,K) 96 : A(I,1)=A(I,1)+CO 97 : A(I,2)=A(I,2)+CC 98 : A(I,3)=A(I,3)+CO*Q(2,K) 99 : A(I,6)=A(I,6)+CC*Q(2,K) 100 : A(I,9)=A(I,9)-CO 101 : A(I,14)=A(I,14)-CO*Q(2,K) 102 : A(I,15)=A(I,15)-CC 103 : 32 A(I,20)=A(I,20)-CC*Q(2,K) 104 : C 105 : DO 15 I=1,3 106 : A(I,1)=A(I,1)-XX(I,IG) 107 : 15 X(I)=0.D0 108 : C 109 : ITER=0 110 : ICRI=0 111 : 100 CONTINUE 112 : C --- LES DERIVEES DES FONCTIONS DE TRANSFORMATION 113 : X11=X(1)*X(1) 114 : X22=X(2)*X(2) 115 : X33=X(3)*X(3) 116 : X12=X(1)*X(2) 117 : X23=X(2)*X(3) 118 : X31=X(3)*X(1) 119 : X123=X(1)*X(2)*X(3) 120 : DO 20 I=1,3 121 : AJ(I,1)=A(I,2)+TWO*A(I,5)*X(1)+A(I,6)*X(2)+A(I,10)*X(3)+ 122 : & TWO*A(I,11)*X12+A(I,12)*X22+A(I,15)*X33+TWO*A(I,16)*X31+ 123 : & A(I,17)*X23+TWO*A(I,18)*X123+A(I,19)*X22*X(3)+A(I,20)*X(2)*X33 124 : C 125 : AJ(I,2)=A(I,3)+A(I,6)*X(1)+TWO*A(I,7)*X(2)+A(I,8)*X(3)+ 126 : & A(I,11)*X11+TWO*A(I,12)*X12+TWO*A(I,13)*X23+A(I,14)*X33+ 127 : & A(I,17)*X31+A(I,18)*X11*X(3)+TWO*A(I,19)*X123+A(I,20)*X(1)*X33 128 : C 129 : 20 AJ(I,3)=A(I,4)+A(I,8)*X(2)+TWO*A(I,9)*X(3)+A(I,10)*X(1)+ 130 : & A(I,13)*X22+TWO*A(I,14)*X23+TWO*A(I,15)*X31+A(I,16)*X11+ 131 : & A(I,17)*X12+A(I,18)*X11*X(2)+A(I,19)*X(1)*X22+TWO*A(I,20)*X123 132 : C --- Calcul de l'inverse de la matrice AJ et son determinant :BJ=inverse(AJ) 132 : CALL INVERT(AJ,BJ,DET,ICRI) 133 : DO 40 I=1,3 134 : 40 F(I)=A(I,1)+A(I,2)*X(1)+A(I,3)*X(2)+A(I,4)*X(3)+A(I,5)*X11 135 : & +A(I,6)*X12+A(I,7)*X22+A(I,8)*X23+A(I,9)*X33+A(I,10)*X31 136 : & +A(I,11)*X11*X(2)+A(I,12)*X(1)*X22+A(I,13)*X22*X(3) 137 : & +A(I,14)*X(2)*X33+A(I,15)*X(1)*X33+A(I,16)*X11*X(3) 138 : & +(A(I,17)+A(I,18)*X(1)+A(I,19)*X(2)+A(I,20)*X(3))*X123 139 : DD=SQRT(F(1)*F(1)+F(2)*F(2)+F(3)*F(3)) 140 : IF(DD.GT.EPSI) THEN 141 : ITER=ITER+1 142 : DO 50 I=1,3 143 : 50 X(I)=X(I)-BJ(I,1)*F(1)-BJ(I,2)*F(2)-BJ(I,3)*F(3) 144 : GOTO 100 145 : END IF 146 : DO 60 I=1,3 147 : 60 QSI(I,IG)=X(I) 148 : C --- Calcul des fonctions d'interpolation leur derivees 149 : C --- Fonctions d'interpolation Ni 150 : DO 310 I=1,4 151 : I1=2*I-1 152 : I2=I1+12 153 : S1=1.D0+X(1)*Q(1,I) 154 : S2=1.D0+X(2)*Q(2,I) 155 : S3=1.D0+X(3)*Q(3,I) 156 : S4=-1.D0+X(1)*Q(1,I)+X(2)*Q(2,I)+X(3)*Q(3,I) 157 : J=I+4 158 : T1=1.D0+X(1)*Q(1,J) 159 : T2=1.D0+X(2)*Q(2,J) 160 : T3=1.D0+X(3)*Q(3,J) 161 : T4=-1.D0+X(1)*Q(1,J)+X(2)*Q(2,J)+X(3)*Q(3,J) 162 : C --- NOEUDS 1 3 5 7 13 15 17 19 163 : SHPC(1,I1,IG)=0.125D0*S1*S2*S3*(-UN+S4) 164 : SHPC(1,I2,IG)=0.125D0*T1*T2*T3*(-UN+T4) 165 : IF(ICM.NE.0) THEN 166 : SHPC(2,I1,IG)=0.125D0*Q(1,I)*S2*S3*(S4+X(1)*Q(1,I)) 167 : SHPC(2,I2,IG)=0.125D0*Q(1,J)*T2*T3*(T4+X(1)*Q(1,J)) 168 : SHPC(3,I1,IG)=0.125D0*Q(2,I)*S1*S3*(S4+X(2)*Q(2,I)) 169 : SHPC(3,I2,IG)=0.125D0*Q(2,J)*T1*T3*(T4+X(2)*Q(2,J)) 170 : SHPC(4,I1,IG)=0.125D0*Q(3,I)*S1*S2*(S4+X(3)*Q(3,I)) 171 : SHPC(4,I2,IG)=0.125D0*Q(3,J)*T1*T2*(T4+X(3)*Q(3,J)) 172 : END IF 173 : C 174 : C --- NOEUDS 2 6 14 18 175 : S0=1.D0-X(1)*X(1) 176 : I3=(I-1)*4+2 177 : IF(I.GE.3) I3=I*4+2 178 : J=I+8 179 : S1=1.D0+X(2)*Q(2,J) 180 : S2=1.D0+X(3)*Q(3,J) 181 : SHPC(1,I3,IG)=0.25D0*S0*S1*S2 182 : IF(ICM.NE.0) THEN 183 : SHPC(2,I3,IG)=-0.5D0*X(1)*S1*S2 184 : SHPC(3,I3,IG)=0.25D0*Q(2,J)*S0*S2 185 : SHPC(4,I3,IG)=0.25D0*Q(3,J)*S0*S1 186 : END IF 187 : C 188 : C --- NOEUDS 4 8 16 20 189 : I4=I*4 190 : IF(I.GE.3) I4=I*4+4 191 : J=I+12 192 : S0=1.D0-X(2)*X(2) 193 : S1=1.D0+X(1)*Q(1,J) 194 : S2=1.D0+X(3)*Q(3,J) 195 : SHPC(1,I4,IG)=0.25D0*S0*S1*S2 196 : IF(ICM.NE.0) THEN 197 : SHPC(2,I4,IG)=0.25D0*Q(1,J)*S0*S2 198 : SHPC(3,I4,IG)=-0.5D0*X(2)*S1*S2 199 : SHPC(4,I4,IG)=0.25D0*Q(3,J)*S0*S1 200 : END IF 201 : C 202 : C --- NOEUDS 9 10 11 12 203 : I5=I+8 204 : J=I+16 205 : S0=1.D0-X(3)*X(3) 206 : S1=1.D0+X(1)*Q(1,J) 207 : S2=1.D0+X(2)*Q(2,J) 208 : SHPC(1,I5,IG)=0.25D0*S0*S1*S2 209 : IF(ICM.NE.0) THEN 210 : SHPC(2,I5,IG)=0.25D0*Q(1,J)*S0*S2 211 : SHPC(3,I5,IG)=0.25D0*Q(2,J)*S0*S1 212 : SHPC(4,I5,IG)=-0.5D0*X(3)*S1*S2 213 : END IF 214 : 310 CONTINUE 214 : CALL PROSGA(ICM,XE,XEC,XX,X,SHPC,SGA,QSI,IG,20) 215 : END 216 : 217 : 218 : 219 : 220 :
© Cast3M 2003 - All rights reserved.
Disclaimer