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

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