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