Télécharger la plaquette Cast3M


fr en

Télécharger qsijs2.eso

Retour à la liste

Numérotation des lignes : 

qsijs2

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

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