C QSIJS2 SOURCE CB215821 19/03/18 21:16:10 10160 SUBROUTINE QSIJS2(XEL,MELE,NBMA1,IDI,XPI,SHP,QSI,IROUT) C======================================================================C C C C INPUT: C C XEL(3,*) : LES COORDONNEES DE L'ELEMENT DANS LE REPERE LOCAL C C MELE : NUMERO DE L ELEMENT C C NBMA1 NOMBRE DE NOEUDS DE L ELEMENT C C XPI(3): COORDONNEES DU POINT DANS LE SYSTEME GLOBAL C C IDI 2 POUR 2D ET ELEMENTS SURFACIQUES 3 POUR ELEM. VOLUMIQUES C C C C INPUT/OUTPUT: C C QSI(3): COORDONNEES DU POINT DANS L ELEMENT DE REFERENCE C C C C OUTPUT: C C C C SHP FONCTIONS DE FORMES ET DERIVEES C C IROUT 0 SI SUCCES DE L OPERATION 1 SI ECHEC C C Sortie convergence a 1.E-5 sur élement de référence C C J.S. FLEURET / TRANSFORME JM BAZE : POUR TOUS ELEMENTS(EN PRINCIPE) C C======================================================================C C Gounand 13/07/2007 : la différence par rapport à qsijs est qu'on C fournit une estimation initiale du tableau QSI IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC CCGEOME DIMENSION XEL(3,*),QSI(3),XPI(3),DQS(3) DIMENSION XSU(20),YSU(20),ZSU(20) DIMENSION SHP(6,*),D(3,3),DINV(3,3),V(3),XU(3),F(3) EXTERNAL SHAPE DATA UN,XZER,EPSI/1.D0,0.D0,1.D-7/ C- Quelques initialisations IROUT=0 Iter=0 IAcc=0 IA1=4 * QSI(1)=XZER * QSI(2)=XZER * QSI(3)=XZER IF (IDI.EQ.1) GOTO 800 DXDQSI=XZER DXDETA=XZER DYDQSI=XZER DYDETA=XZER C-------------------- DEBUT DES ITERATIONS ------------------------ 600 CONTINUE ITER = ITER + 1 IACC = IACC + 1 IF (ITER.GT.100) THEN IROUT = ITER RETURN ENDIF * WRITE(IOIMP,*) ' qsijs2 : QSI1=',QSI(1),' QSI2=',QSI(2) CALL SHAPE(QSI(1),QSI(2),QSI(3),MELE,SHP,IRET) if (iret.eq.0) then MOTERR(1:4)=NOMS(MELE) CALL ERREUR(68) RETURN endif C PROJECTION DO 501 IK = 1,IDI XU(IK)= XZER DO 501 JK = 1,NBMA1 XU(IK) = XU(IK) + SHP(1,JK)*(XEL(IK,JK)) 501 CONTINUE C DD = XZER DO 502 IK = 1,IDI F(IK) = XU(IK) - XPI(IK) 502 CONTINUE C C CAS 2 ELEMENTS SURFACIQUES C IF (IDI.EQ.2) THEN C IF ((MELE.LT.14).OR.(MELE.GT.26)) THEN DO 100 I=1,NBMA1 DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I) DXDETA=DXDETA+SHP(3,I)*XEL(1,I) DYDQSI=DYDQSI+SHP(2,I)*XEL(2,I) DYDETA=DYDETA+SHP(3,I)*XEL(2,I) 100 CONTINUE DJAC=DXDQSI*DYDETA-DXDETA*DYDQSI XXXX = DJAC IF(DJAC.NE.XZER) XXXX=UN/DJAC DQSIDX= DYDETA*XXXX DQSIDY=-DXDETA*XXXX DETADX=-DYDQSI*XXXX DETADY= DXDQSI*XXXX DINV(1,1) = DQSIDX DINV(2,1) = DQSIDY DINV(1,2) = DETADX DINV(2,2) = DETADY ELSE C C ELEMENTS VOLUMIQUES C DO 310 I=1,3 DO 310 J=1,3 D(I,J)=XZER 310 CONTINUE DO 300 I=1,NBMA1 C D(1,1)=D(1,1)+SHP(2,I)*XEL(1,I) D(2,1)=D(2,1)+SHP(3,I)*XEL(1,I) D(3,1)=D(3,1)+SHP(4,I)*XEL(1,I) C D(1,2)=D(1,2)+SHP(2,I)*XEL(2,I) D(2,2)=D(2,2)+SHP(3,I)*XEL(2,I) D(3,2)=D(3,2)+SHP(4,I)*XEL(2,I) C D(1,3)=D(1,3)+SHP(2,I)*XEL(3,I) D(2,3)=D(2,3)+SHP(3,I)*XEL(3,I) D(3,3)=D(3,3)+SHP(4,I)*XEL(3,I) 300 CONTINUE C INVERSION DE LA MATRICE D(3,3) DINV(1,1)= D(2,2)*D(3,3)-D(2,3)*D(3,2) DINV(2,2)= D(1,1)*D(3,3)-D(1,3)*D(3,1) DINV(3,3)= D(1,1)*D(2,2)-D(1,2)*D(2,1) DINV(1,2)=-D(1,2)*D(3,3)+D(3,2)*D(1,3) DINV(2,1)=-D(2,1)*D(3,3)+D(2,3)*D(3,1) DINV(1,3)= D(1,2)*D(2,3)-D(2,2)*D(1,3) DINV(3,1)= D(2,1)*D(3,2)-D(2,2)*D(3,1) DINV(2,3)=-D(1,1)*D(2,3)+D(2,1)*D(1,3) DINV(3,2)=-D(1,1)*D(3,2)+D(1,2)*D(3,1) DJAC=D(1,1)*DINV(1,1)+D(2,1)*DINV(1,2)+D(3,1)*DINV(1,3) XXXX = DJAC IF(DJAC.NE.XZER) XXXX=UN/DJAC DO 350 IA=1,3 DO 350 IB=1,3 DINV(IA,IB)=DINV(IA,IB)*XXXX 350 CONTINUE ENDIF C DDQ =XZER DO 505 I=1,IDI DQS(I)=XZER DO 506 J=1,IDI DQS(I) =DQS(I)+ DINV(J,I)*F(J) 506 CONTINUE DDQ = DDQ + DQS(I)*DQS(I) 505 QSI(I) = QSI(I)-DQS(I) DDQ = SQRT(DDQ) IF(DDQ.LT.EPSI) GOTO 699 XSU(IACC) = QSI(1) YSU(IACC) = QSI(2) IF(IDI.EQ.3) ZSU(IACC) = QSI(3) C IF(IACC.GT.IA1 .AND. IACC.GT.2) THEN C acceleration dqs1 = 2.D0*XSU(IACC-1)-QSI(1)-XSU(IACC-2) if (abs(dqs1).GT.1.D-10) then BX = (QSI(1) - XSU(IACC-1))*2.D0/dqs1 QSI(1) =(1.D0+BX)*XSU(IACC-1) - BX*XSU(IACC-2) end if dqs2 = 2.D0*YSU(IACC-1)-QSI(2)-YSU(IACC-2) if (abs(dqs2).GT.1.D-10) then BY = (QSI(2) - YSU(IACC-1))*2.D0/dqs2 QSI(2) =(1.D0+BY)*YSU(IACC-1) - BY*YSU(IACC-2) endif XSU(IACC) = QSI(1) YSU(IACC) = QSI(2) IF(IDI.EQ.3) THEN dqs3 = 2.D0*ZSU(IACC-1)-QSI(3)-ZSU(IACC-2) if (abs(dqs3).GT.1.D-10) then BZ = (QSI(3) - ZSU(IACC-1))*2.D0/dqs3 QSI(3) =(1.D0+BZ)*ZSU(IACC-1) - BZ*ZSU(IACC-2) endif ZSU(IACC) = QSI(3) ENDIF IA1 = IA1+2 ENDIF C au cas ou il faudrait plus de 20 acceleration if ( iacc.EQ.20) then do 400 i = 1,2 xsu(i) = xsu(iacc-2+i) ysu(i) = ysu(iacc-2+i) zsu(i) = zsu(iacc-2+i) 400 continue ia1 = 1 iacc =0 endif C C-------- ITERATIONS ---------------------- GOTO 600 * ENDIF C========== Cas IDI = 1 ========== C Methode de Newton-Raphson pour resoudre X(QSI)-XPI(1)=0 800 Iter=Iter+1 IAcc=IAcc+1 IF (Iter.GT.100) THEN IROUT=Iter RETURN ENDIF CALL SHAPE(QSI(1),QSI(2),QSI(3),MELE,SHP,IRET) if (iret.eq.0) then MOTERR(1:4)=NOMS(MELE) CALL ERREUR(68) RETURN endif XU1=XZER dXdQSI=XZER DO i=1,NBMA1 XU1=XU1+SHP(1,i)*XEL(1,i) dXdQSI=dXdQSI+SHP(2,i)*XEL(1,i) ENDDO IF (dXdQSI.NE.XZER) THEN F1=XU1-XPI(1) dQS1=F1/dXdQSI QSI(1)=QSI(1)-dQS1 ELSE dQS1=XZER ENDIF C Test de convergence IF (ABS(dQS1).LT.EPSI) RETURN C Acceleration IF (IAcc.GT.IA1) THEN dQS1=2.*XSU(IAcc-1)-QSI(1)-XSU(IAcc-2) IF (ABS(dQS1).GT.1.D-10) THEN BX=2.*(QSI(1)-XSU(IAcc-1))/dQS1 QSI(1)=(UN+BX)*XSU(IAcc-1)-BX*XSU(IAcc-2) ENDIF IA1=IA1+2 ENDIF XSU(IAcc)=QSI(1) C Dans le cas ou il faudrait plus de 20 accelerations IF (IAcc.EQ.20) THEN XSU(1)=XSU(IAcc-1) XSU(2)=XSU(IAcc) IA1=1 IAcc=0 ENDIF GOTO 800 699 CONTINUE END