Télécharger la plaquette Cast3M


fr en

Télécharger qsijs.eso

Retour à la liste

Numérotation des lignes : 

qsijs

  1. C QSIJS     SOURCE    CB215821  19/09/30    21:15:20     10324          
  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.  
  22.  
  23. -INC PPARAM
  24. -INC CCOPTIO
  25. -INC CCGEOME
  26.  
  27.       DIMENSION XEL(3,*),QSI(3),XPI(3)
  28.       DIMENSION XSU(20),YSU(20),ZSU(20)
  29.       DIMENSION SHP(6,*),D(3,3),DINV(3,3),V(3),XU(3),F(3)
  30.       EXTERNAL SHAPE
  31.       DATA UN,XZER,EPSI/1.D0,0.D0,1.D-7/
  32.  
  33. C- Quelques initialisations
  34.       IROUT=0
  35.       Iter=0
  36.       IAcc=0
  37.       IA1=4
  38.  
  39.       QSI(1)=XZER
  40.       QSI(2)=XZER
  41.       QSI(3)=XZER
  42.  
  43.       IF (IDI.EQ.1) GOTO 800
  44.  
  45.       DXDQSI=XZER
  46.       DXDETA=XZER
  47.       DYDQSI=XZER
  48.       DYDETA=XZER
  49.  
  50. C-------------------- DEBUT DES ITERATIONS ------------------------
  51. 600   CONTINUE
  52.       ITER = ITER + 1
  53.       IACC = IACC + 1
  54.       IF (ITER.GT.100) THEN
  55.       IROUT = ITER
  56.       RETURN
  57.       ENDIF
  58. *      WRITE(IOIMP,*) ' qsijs : 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)
  62.          CALL ERREUR(68)
  63.          RETURN
  64.       endif
  65.  
  66.  
  67. C           PROJECTION
  68.       CALL ZERO(XU,IDI,1)
  69.       DO 501 JK = 1,NBMA1
  70.         SHPK=SHP(1,JK)
  71.         DO 502 IK = 1,IDI
  72.           XU(IK) = XU(IK) + SHPK*XEL(IK,JK)
  73. 502     CONTINUE
  74. 501   CONTINUE
  75.  
  76.       DO 503 IK = 1,IDI
  77.         F(IK) = XU(IK) - XPI(IK)
  78. 503   CONTINUE
  79. C
  80. C     CAS 2 ELEMENTS SURFACIQUES
  81. C
  82.       IF (IDI.EQ.2) THEN
  83. C      IF ((MELE.LT.14).OR.(MELE.GT.26)) THEN
  84.       DO 100 I=1,NBMA1
  85.       DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I)
  86.       DXDETA=DXDETA+SHP(3,I)*XEL(1,I)
  87.       DYDQSI=DYDQSI+SHP(2,I)*XEL(2,I)
  88.       DYDETA=DYDETA+SHP(3,I)*XEL(2,I)
  89.  100  CONTINUE
  90.       DJAC=DXDQSI*DYDETA-DXDETA*DYDQSI
  91.       XXXX = DJAC
  92.       IF(DJAC.NE.XZER) XXXX=UN/DJAC
  93.       DQSIDX= DYDETA*XXXX
  94.       DQSIDY=-DXDETA*XXXX
  95.       DETADX=-DYDQSI*XXXX
  96.       DETADY= DXDQSI*XXXX
  97.       DINV(1,1) = DQSIDX
  98.       DINV(2,1) = DQSIDY
  99.       DINV(1,2) = DETADX
  100.       DINV(2,2) = DETADY
  101.  
  102.       ELSE
  103. C
  104. C     ELEMENTS VOLUMIQUES
  105. C
  106.       CALL ZERO(D,3,3)
  107.       DO 300 I=1,NBMA1
  108.         SHP2=SHP(2,I)
  109.         SHP3=SHP(3,I)
  110.         SHP4=SHP(4,I)
  111.  
  112.         XELi=XEL(1,I)
  113.         D(1,1)=D(1,1)+SHP2*XELi
  114.         D(2,1)=D(2,1)+SHP3*XELi
  115.         D(3,1)=D(3,1)+SHP4*XELi
  116. C
  117.         XELi=XEL(2,I)
  118.         D(1,2)=D(1,2)+SHP2*XELi
  119.         D(2,2)=D(2,2)+SHP3*XELi
  120.         D(3,2)=D(3,2)+SHP4*XELi
  121. C
  122.         XELi=XEL(3,I)
  123.         D(1,3)=D(1,3)+SHP2*XELi
  124.         D(2,3)=D(2,3)+SHP3*XELi
  125.         D(3,3)=D(3,3)+SHP4*XELi
  126.  300  CONTINUE
  127.  
  128.  
  129. C     INVERSION DE LA MATRICE D(3,3)
  130.       DINV(1,1)= D(2,2)*D(3,3)-D(2,3)*D(3,2)
  131.       DINV(2,2)= D(1,1)*D(3,3)-D(1,3)*D(3,1)
  132.       DINV(3,3)= D(1,1)*D(2,2)-D(1,2)*D(2,1)
  133.       DINV(1,2)=-D(1,2)*D(3,3)+D(3,2)*D(1,3)
  134.       DINV(2,1)=-D(2,1)*D(3,3)+D(2,3)*D(3,1)
  135.       DINV(1,3)= D(1,2)*D(2,3)-D(2,2)*D(1,3)
  136.       DINV(3,1)= D(2,1)*D(3,2)-D(2,2)*D(3,1)
  137.       DINV(2,3)=-D(1,1)*D(2,3)+D(2,1)*D(1,3)
  138.       DINV(3,2)=-D(1,1)*D(3,2)+D(1,2)*D(3,1)
  139.       DJAC=D(1,1)*DINV(1,1)+D(2,1)*DINV(1,2)+D(3,1)*DINV(1,3)
  140.       XXXX = DJAC
  141.       IF(DJAC.NE.XZER) XXXX=UN/DJAC
  142.       DO 350 IA=1,3
  143.       DO 350 IB=1,3
  144.       DINV(IB,IA)=DINV(IB,IA)*XXXX
  145.   350 CONTINUE
  146.       ENDIF
  147. C
  148.        DDQ =XZER
  149.        DO 505 I=1,IDI
  150.          XFLO  =XZER
  151.          DO 506 J=1,IDI
  152.            XFLO =XFLO + DINV(J,I)*F(J)
  153. 506      CONTINUE
  154.          DDQ    = DDQ + (XFLO ** 2)
  155.          QSI(I) = QSI(I)-XFLO
  156. 505    CONTINUE
  157.        DDQ = SQRT(DDQ)
  158.        IF(DDQ.LT.EPSI) GOTO 699
  159.        XSU(IACC) = QSI(1)
  160.        YSU(IACC) = QSI(2)
  161.        IF(IDI.EQ.3) ZSU(IACC) = QSI(3)
  162. C
  163.          IF(IACC.GT.IA1 .AND. IACC.GT.2) THEN
  164. C          acceleration
  165.        dqs1 = 2.D0*XSU(IACC-1)-QSI(1)-XSU(IACC-2)
  166.        if (abs(dqs1).GT.1.D-10) then
  167.        BX = (QSI(1) - XSU(IACC-1))*2.D0/dqs1
  168.        QSI(1) =(1.D0+BX)*XSU(IACC-1) - BX*XSU(IACC-2)
  169.        end if
  170.        dqs2 = 2.D0*YSU(IACC-1)-QSI(2)-YSU(IACC-2)
  171.        if (abs(dqs2).GT.1.D-10) then
  172.        BY = (QSI(2) - YSU(IACC-1))*2.D0/dqs2
  173.        QSI(2) =(1.D0+BY)*YSU(IACC-1) - BY*YSU(IACC-2)
  174.        endif
  175.        XSU(IACC) = QSI(1)
  176.        YSU(IACC) = QSI(2)
  177.          IF(IDI.EQ.3) THEN
  178.        dqs3 = 2.D0*ZSU(IACC-1)-QSI(3)-ZSU(IACC-2)
  179.          if (abs(dqs3).GT.1.D-10) then
  180.        BZ = (QSI(3) - ZSU(IACC-1))*2.D0/dqs3
  181.        QSI(3) =(1.D0+BZ)*ZSU(IACC-1) - BZ*ZSU(IACC-2)
  182.          endif
  183.        ZSU(IACC) = QSI(3)
  184.          ENDIF
  185.          IA1 = IA1+2
  186.        ENDIF
  187. C   au cas  ou il faudrait plus de 20 acceleration
  188.        if ( iacc.EQ.20) then
  189.         do 400 i = 1,2
  190.        xsu(i) = xsu(iacc-2+i)
  191.        ysu(i) = ysu(iacc-2+i)
  192.        zsu(i) = zsu(iacc-2+i)
  193. 400     continue
  194.        ia1 = 1
  195.        iacc =0
  196.        endif
  197.  
  198. C
  199. C-------- ITERATIONS ----------------------
  200.  
  201.        GOTO 600
  202. *      ENDIF
  203.  
  204. C==========  Cas IDI = 1  ==========
  205. C Methode de Newton-Raphson pour resoudre X(QSI)-XPI(1)=0
  206.  800  Iter=Iter+1
  207.       IAcc=IAcc+1
  208.       IF (Iter.GT.100) THEN
  209.         IROUT=Iter
  210.         RETURN
  211.       ENDIF
  212.       CALL SHAPE(QSI(1),QSI(2),QSI(3),MELE,SHP,IRET)
  213.       if (iret.eq.0) then
  214.          MOTERR(1:4)=NOMS(MELE)
  215.          CALL ERREUR(68)
  216.          RETURN
  217.       endif
  218.       XU1=XZER
  219.       dXdQSI=XZER
  220.       DO i=1,NBMA1
  221.         XU1=XU1+SHP(1,i)*XEL(1,i)
  222.         dXdQSI=dXdQSI+SHP(2,i)*XEL(1,i)
  223.       ENDDO
  224.       IF (dXdQSI.NE.XZER) THEN
  225.         F1=XU1-XPI(1)
  226.         dQS1=F1/dXdQSI
  227.         QSI(1)=QSI(1)-dQS1
  228.       ELSE
  229.         dQS1=XZER
  230.       ENDIF
  231. C Test de convergence
  232.       IF (ABS(dQS1).LT.EPSI) RETURN
  233. C Acceleration
  234.       IF (IAcc.GT.IA1) THEN
  235.         dQS1=2.*XSU(IAcc-1)-QSI(1)-XSU(IAcc-2)
  236.         IF (ABS(dQS1).GT.1.D-10) THEN
  237.           BX=2.*(QSI(1)-XSU(IAcc-1))/dQS1
  238.           QSI(1)=(UN+BX)*XSU(IAcc-1)-BX*XSU(IAcc-2)
  239.         ENDIF
  240.         IA1=IA1+2
  241.       ENDIF
  242.       XSU(IAcc)=QSI(1)
  243. C Dans le cas ou il faudrait plus de 20 accelerations
  244.       IF (IAcc.EQ.20) THEN
  245.         XSU(1)=XSU(IAcc-1)
  246.         XSU(2)=XSU(IAcc)
  247.         IA1=1
  248.         IAcc=0
  249.       ENDIF
  250.       GOTO 800
  251.  
  252.  699  CONTINUE
  253.       END
  254.  
  255.  
  256.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales