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

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