1 : C QSIJS SOURCE CHAT 07/10/29 21:27:08 5940 2 : SUBROUTINE QSIJS(XEL,MELE,NBMA1,IDI,XPI,SHP,QSI,IROUT) 3 : C======================================================================C 4 : C C 5 : C INPUT: C 6 : C XEL(3,*) : LES COORDONNEES DE L'ELEMENT DANS LE REPERE LOCAL C C 7 : C MELE : NUMERO DE L ELEMENT C 8 : C NBMA1 NOMBRE DE NOEUDS DE L ELEMENT C 9 : C XPI(3): COORDONNEES DU POINT DANS LE SYSTEME GLOBAL C 10 : C IDI 2 POUR 2D ET ELEMENTS SURFACIQUES 3 POUR ELEM. VOLUMIQUES C 11 : C OUTPUT: C 12 : C C 13 : C QSI(3): COORDONNEES DU POINT DANS L ELEMENT DE REFERENCE C 14 : C SHP FONCTIONS DE FORMES ET DERIVEES C 15 : C IROUT 0 SI SUCCES DE L OPERATION 1 SI ECHEC C C 16 : C Sortie convergence a 1.E-5 sur élement de référence C C 17 : C J.S. FLEURET / TRANSFORME JM BAZE : POUR TOUS ELEMENTS(EN PRINCIPE) C C 18 : C======================================================================C 19 : IMPLICIT INTEGER(I-N) 20 : IMPLICIT REAL*8 (A-H,O-Z) 21 : -INC CCOPTIO 22 : -INC CCGEOME 23 : DIMENSION XEL(3,*),QSI(3),XPI(3),DQS(3) 24 : DIMENSION XSU(20),YSU(20),ZSU(20) 25 : DIMENSION SHP(6,*),D(3,3),DINV(3,3),V(3),XU(3),F(3) 26 : EXTERNAL SHAPE 27 : DATA UN,XZER,EPSI/1.D0,0.D0,1.D-7/ 28 : 29 : C- Quelques initialisations 30 : IROUT=0 31 : Iter=0 32 : IAcc=0 33 : IA1=4 34 : 35 : QSI(1)=XZER 36 : QSI(2)=XZER 37 : QSI(3)=XZER 38 : 39 : IF (IDI.EQ.1) GOTO 800 40 : 41 : DXDQSI=XZER 42 : DXDETA=XZER 43 : DYDQSI=XZER 44 : DYDETA=XZER 45 : 46 : C-------------------- DEBUT DES ITERATIONS ------------------------ 47 : 600 CONTINUE 48 : ITER = ITER + 1 49 : IACC = IACC + 1 50 : IF (ITER.GT.100) THEN 51 : IROUT = ITER 52 : RETURN 53 : ENDIF 54 : * WRITE(IOIMP,*) ' qsijs : QSI1=',QSI(1),' QSI2=',QSI(2) 54 : CALL SHAPE(QSI(1),QSI(2),QSI(3),MELE,SHP,IRET) 55 : if (iret.eq.0) then 56 : MOTERR(1:4)=NOMS(MELE) 56 : CALL ERREUR(68) 57 : RETURN 58 : endif 59 : 60 : 61 : C PROJECTION 62 : DO 501 IK = 1,IDI 63 : XU(IK)= XZER 64 : DO 501 JK = 1,NBMA1 65 : XU(IK) = XU(IK) + SHP(1,JK)*(XEL(IK,JK)) 66 : 501 CONTINUE 67 : 68 : C DD = XZER 69 : DO 502 IK = 1,IDI 70 : F(IK) = XU(IK) - XPI(IK) 71 : 502 CONTINUE 72 : C 73 : C CAS 2 ELEMENTS SURFACIQUES 74 : C 75 : IF (IDI.EQ.2) THEN 76 : C IF ((MELE.LT.14).OR.(MELE.GT.26)) THEN 77 : DO 100 I=1,NBMA1 78 : DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I) 79 : DXDETA=DXDETA+SHP(3,I)*XEL(1,I) 80 : DYDQSI=DYDQSI+SHP(2,I)*XEL(2,I) 81 : DYDETA=DYDETA+SHP(3,I)*XEL(2,I) 82 : 100 CONTINUE 83 : DJAC=DXDQSI*DYDETA-DXDETA*DYDQSI 84 : XXXX = DJAC 85 : IF(DJAC.NE.XZER) XXXX=UN/DJAC 86 : DQSIDX= DYDETA*XXXX 87 : DQSIDY=-DXDETA*XXXX 88 : DETADX=-DYDQSI*XXXX 89 : DETADY= DXDQSI*XXXX 90 : DINV(1,1) = DQSIDX 91 : DINV(2,1) = DQSIDY 92 : DINV(1,2) = DETADX 93 : DINV(2,2) = DETADY 94 : 95 : ELSE 96 : C 97 : C ELEMENTS VOLUMIQUES 98 : C 99 : 100 : DO 310 I=1,3 101 : DO 310 J=1,3 102 : D(I,J)=XZER 103 : 310 CONTINUE 104 : DO 300 I=1,NBMA1 105 : C 106 : D(1,1)=D(1,1)+SHP(2,I)*XEL(1,I) 107 : D(2,1)=D(2,1)+SHP(3,I)*XEL(1,I) 108 : D(3,1)=D(3,1)+SHP(4,I)*XEL(1,I) 109 : C 110 : D(1,2)=D(1,2)+SHP(2,I)*XEL(2,I) 111 : D(2,2)=D(2,2)+SHP(3,I)*XEL(2,I) 112 : D(3,2)=D(3,2)+SHP(4,I)*XEL(2,I) 113 : C 114 : D(1,3)=D(1,3)+SHP(2,I)*XEL(3,I) 115 : D(2,3)=D(2,3)+SHP(3,I)*XEL(3,I) 116 : D(3,3)=D(3,3)+SHP(4,I)*XEL(3,I) 117 : 300 CONTINUE 118 : C INVERSION DE LA MATRICE D(3,3) 119 : DINV(1,1)= D(2,2)*D(3,3)-D(2,3)*D(3,2) 120 : DINV(2,2)= D(1,1)*D(3,3)-D(1,3)*D(3,1) 121 : DINV(3,3)= D(1,1)*D(2,2)-D(1,2)*D(2,1) 122 : DINV(1,2)=-D(1,2)*D(3,3)+D(3,2)*D(1,3) 123 : DINV(2,1)=-D(2,1)*D(3,3)+D(2,3)*D(3,1) 124 : DINV(1,3)= D(1,2)*D(2,3)-D(2,2)*D(1,3) 125 : DINV(3,1)= D(2,1)*D(3,2)-D(2,2)*D(3,1) 126 : DINV(2,3)=-D(1,1)*D(2,3)+D(2,1)*D(1,3) 127 : DINV(3,2)=-D(1,1)*D(3,2)+D(1,2)*D(3,1) 128 : DJAC=D(1,1)*DINV(1,1)+D(2,1)*DINV(1,2)+D(3,1)*DINV(1,3) 129 : XXXX = DJAC 130 : IF(DJAC.NE.XZER) XXXX=UN/DJAC 131 : DO 350 IA=1,3 132 : DO 350 IB=1,3 133 : DINV(IA,IB)=DINV(IA,IB)*XXXX 134 : 350 CONTINUE 135 : ENDIF 136 : C 137 : DDQ =XZER 138 : DO 505 I=1,IDI 139 : DQS(I)=XZER 140 : DO 506 J=1,IDI 141 : DQS(I) =DQS(I)+ DINV(J,I)*F(J) 142 : 506 CONTINUE 143 : DDQ = DDQ + DQS(I)*DQS(I) 144 : 505 QSI(I) = QSI(I)-DQS(I) 145 : DDQ = SQRT(DDQ) 146 : IF(DDQ.LT.EPSI) GOTO 699 147 : XSU(IACC) = QSI(1) 148 : YSU(IACC) = QSI(2) 149 : IF(IDI.EQ.3) ZSU(IACC) = QSI(3) 150 : C 151 : IF(IACC.GT.IA1) THEN 152 : C acceleration 153 : dqs1 = 2.D0*XSU(IACC-1)-QSI(1)-XSU(IACC-2) 154 : if (abs(dqs1).GT.1.D-10) then 155 : BX = (QSI(1) - XSU(IACC-1))*2.D0/dqs1 156 : QSI(1) =(1.D0+BX)*XSU(IACC-1) - BX*XSU(IACC-2) 157 : end if 158 : dqs2 = 2.D0*YSU(IACC-1)-QSI(2)-YSU(IACC-2) 159 : if (abs(dqs2).GT.1.D-10) then 160 : BY = (QSI(2) - YSU(IACC-1))*2.D0/dqs2 161 : QSI(2) =(1.D0+BY)*YSU(IACC-1) - BY*YSU(IACC-2) 162 : endif 163 : XSU(IACC) = QSI(1) 164 : YSU(IACC) = QSI(2) 165 : IF(IDI.EQ.3) THEN 166 : dqs3 = 2.D0*ZSU(IACC-1)-QSI(3)-ZSU(IACC-2) 167 : if (abs(dqs3).GT.1.D-10) then 168 : BZ = (QSI(3) - ZSU(IACC-1))*2.D0/dqs3 169 : QSI(3) =(1.D0+BZ)*ZSU(IACC-1) - BZ*ZSU(IACC-2) 170 : endif 171 : ZSU(IACC) = QSI(3) 172 : ENDIF 173 : IA1 = IA1+2 174 : ENDIF 175 : C au cas ou il faudrait plus de 20 acceleration 176 : if ( iacc.EQ.20) then 177 : do 400 i = 1,2 178 : xsu(i) = xsu(iacc-2+i) 179 : ysu(i) = ysu(iacc-2+i) 180 : zsu(i) = zsu(iacc-2+i) 181 : 400 continue 182 : ia1 = 1 183 : iacc =0 184 : endif 185 : 186 : C 187 : C-------- ITERATIONS ---------------------- 188 : 189 : GOTO 600 190 : * ENDIF 191 : 192 : C========== Cas IDI = 1 ========== 193 : C Methode de Newton-Raphson pour resoudre X(QSI)-XPI(1)=0 194 : 800 Iter=Iter+1 195 : IAcc=IAcc+1 196 : IF (Iter.GT.100) THEN 197 : IROUT=Iter 198 : RETURN 199 : ENDIF 199 : CALL SHAPE(QSI(1),QSI(2),QSI(3),MELE,SHP,IRET) 200 : if (iret.eq.0) then 201 : MOTERR(1:4)=NOMS(MELE) 201 : CALL ERREUR(68) 202 : RETURN 203 : endif 204 : XU1=XZER 205 : dXdQSI=XZER 206 : DO i=1,NBMA1 207 : XU1=XU1+SHP(1,i)*XEL(1,i) 208 : dXdQSI=dXdQSI+SHP(2,i)*XEL(1,i) 209 : ENDDO 210 : IF (dXdQSI.NE.XZER) THEN 211 : F1=XU1-XPI(1) 212 : dQS1=F1/dXdQSI 213 : QSI(1)=QSI(1)-dQS1 214 : ELSE 215 : dQS1=XZER 216 : ENDIF 217 : C Test de convergence 218 : IF (ABS(dQS1).LT.EPSI) RETURN 219 : C Acceleration 220 : IF (IAcc.GT.IA1) THEN 221 : dQS1=2.*XSU(IAcc-1)-QSI(1)-XSU(IAcc-2) 222 : IF (ABS(dQS1).GT.1.D-10) THEN 223 : BX=2.*(QSI(1)-XSU(IAcc-1))/dQS1 224 : QSI(1)=(UN+BX)*XSU(IAcc-1)-BX*XSU(IAcc-2) 225 : ENDIF 226 : IA1=IA1+2 227 : ENDIF 228 : XSU(IAcc)=QSI(1) 229 : C Dans le cas ou il faudrait plus de 20 accelerations 230 : IF (IAcc.EQ.20) THEN 231 : XSU(1)=XSU(IAcc-1) 232 : XSU(2)=XSU(IAcc) 233 : IA1=1 234 : IAcc=0 235 : ENDIF 236 : GOTO 800 237 : 238 : 699 CONTINUE 239 : RETURN 240 : END 241 : 242 : 243 : 244 : 245 :
© Cast3M 2003 - All rights reserved.
Disclaimer