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  25/07/28    21:15:10     12336          
  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     QSI(3)   : COORDONNEES INITIALES DE DEBUT DE RECHERCHE           C
  12. C                                                                      C
  13. C     OUTPUT:                                                          C
  14. C     QSI(3)   : COORDONNEES DU POINT DANS L ELEMENT DE REFERENCE      C
  15. C     SHP      : FONCTIONS DE FORMES ET DERIVEES                       C
  16. C     IROUT   0 SI  SUCCES DE L OPERATION   1  SI ECHEC                C
  17. C        Sortie  convergence a 1.E-7 sur element de reference          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.  
  31.       PARAMETER ( NITMAX = 100 )
  32.       DIMENSION D(3,3),DINV(3,3)
  33. C Pour l'acceleration de convergence (Delta2 Aitken) :
  34.       PARAMETER ( ITAC1 = 5 , ITACN = 2 )
  35. C      DIMENSION XSU(1+ITACN),YSU(1+ITACN),ZSU(1+ITACN)
  36.       DIMENSION XSU(3),YSU(3),ZSU(3)
  37.  
  38.       EXTERNAL SHAPE
  39.       DATA XUN,XZER,EPSI/1.D0,0.D0,1.D-7/
  40.  
  41. C- Quelques initialisations
  42.       IROUT = 0
  43.       ITER  = 0
  44.  
  45. C- Le point de depart (initialisation) est donne avant l'appel.
  46.  
  47.       ITACC = ITAC1
  48.       IACC = 1
  49.       IAC1 = 3
  50.       IAC2 = 2
  51.  
  52.       IF (IDI.EQ.3) GOTO 600
  53.       IF (IDI.EQ.2) GOTO 700
  54.       IF (IDI.EQ.1) GOTO 800
  55.       write(ioimp,*) 'ERREUR INACCEPTABLE DANS QSIJS !'
  56.       call erreur(5)
  57.       irout = 9995
  58.       goto 900
  59.  
  60. C==========  Cas IDI = 3  ==============================================
  61.  600  CONTINUE
  62. C
  63. C     ELEMENTS VOLUMIQUES
  64. C
  65.       EPSIp2 = EPSI*EPSI
  66.  
  67. C-------------------- DEBUT DES ITERATIONS ------------------------
  68.  610  CONTINUE
  69.         ITER = ITER + 1
  70.         IF (ITER.GT.NITMAX) THEN
  71.           IROUT = ITER
  72.           GOTO 900
  73.         ENDIF
  74.  
  75.         CALL SHAPE(QSI(1),QSI(2),QSI(3),MELE,SHP,iret)
  76.         if (iret.eq.0) then
  77.           moterr(1:4) = noms(mele)
  78.           call erreur(68)
  79.           irout = 9999
  80.           goto 900
  81.         endif
  82.  
  83.         F_1 = XZER
  84.         F_2 = XZER
  85.         F_3 = XZER
  86.         DO ik = 1, NBMA1
  87.           F_1 = F_1 + SHP(1,ik) * XEL(1,ik)
  88.           F_2 = F_2 + SHP(1,ik) * XEL(2,ik)
  89.           F_3 = F_3 + SHP(1,ik) * XEL(3,ik)
  90.         ENDDO
  91.         F_1 = F_1 - XPG(1)
  92.         F_2 = F_2 - XPG(2)
  93.         F_3 = F_3 - XPG(3)
  94.  
  95. c*        F_DIST = (F_1 * F_1) + (F_2 * F_2) + (F_3 * F_3)
  96. c*        IF (F_DIST .LT. EPSIp2) GOTO 900
  97.  
  98.         CALL ZERO(D,3,3)
  99.         DO ik = 1, NBMA1
  100.           SHP2=SHP(2,ik)
  101.           SHP3=SHP(3,ik)
  102.           SHP4=SHP(4,ik)
  103.  
  104.           XELi=XEL(1,ik)
  105.           D(1,1)=D(1,1)+SHP2*XELi
  106.           D(2,1)=D(2,1)+SHP3*XELi
  107.           D(3,1)=D(3,1)+SHP4*XELi
  108.  
  109.           XELi=XEL(2,ik)
  110.           D(1,2)=D(1,2)+SHP2*XELi
  111.           D(2,2)=D(2,2)+SHP3*XELi
  112.           D(3,2)=D(3,2)+SHP4*XELi
  113.  
  114.           XELi=XEL(3,ik)
  115.           D(1,3)=D(1,3)+SHP2*XELi
  116.           D(2,3)=D(2,3)+SHP3*XELi
  117.           D(3,3)=D(3,3)+SHP4*XELi
  118.         ENDDO
  119.  
  120. C     INVERSION DE LA MATRICE D(3,3)
  121.         DINV(1,1)= D(2,2)*D(3,3)-D(2,3)*D(3,2)
  122.         DINV(2,2)= D(1,1)*D(3,3)-D(1,3)*D(3,1)
  123.         DINV(3,3)= D(1,1)*D(2,2)-D(1,2)*D(2,1)
  124.         DINV(1,2)=-D(1,2)*D(3,3)+D(3,2)*D(1,3)
  125.         DINV(2,1)=-D(2,1)*D(3,3)+D(2,3)*D(3,1)
  126.         DINV(1,3)= D(1,2)*D(2,3)-D(2,2)*D(1,3)
  127.         DINV(3,1)= D(2,1)*D(3,2)-D(2,2)*D(3,1)
  128.         DINV(2,3)=-D(1,1)*D(2,3)+D(2,1)*D(1,3)
  129.         DINV(3,2)=-D(1,1)*D(3,2)+D(1,2)*D(3,1)
  130.         DJAC=D(1,1)*DINV(1,1)+D(2,1)*DINV(1,2)+D(3,1)*DINV(1,3)
  131.         XX = DJAC
  132.         IF (DJAC.NE.XZER) XX = XUN / DJAC
  133.  
  134.         dqs1 = XX * (DINV(1,1)*F_1 + DINV(2,1)*F_2 + DINV(3,1)*F_3)
  135.         dqs2 = XX * (DINV(1,2)*F_1 + DINV(2,2)*F_2 + DINV(3,2)*F_3)
  136.         dqs3 = XX * (DINV(1,3)*F_1 + DINV(2,3)*F_2 + DINV(3,3)*F_3)
  137.         QSI(1) = QSI(1) - dqs1
  138.         QSI(2) = QSI(2) - dqs2
  139.         QSI(3) = QSI(3) - dqs3
  140.         dDQ = (dqs1 * dqs1) + (dqs2 * dqs2) + (dqs3 * dqs3)
  141.         IF (dDQ .LT. EPSIp2) GOTO 900
  142.  
  143. C Acceleration de convergence :
  144.         IF (ITER.EQ.ITACC) THEN
  145.           dqs1 = 2.D0*XSU(IAC1)-QSI(1)-XSU(IAC2)
  146.           if (abs(dqs1).GT.1.D-10) then
  147.             BX = (QSI(1) - XSU(IAC1))*2.D0/dqs1
  148.             QSI(1) =(1.D0+BX)*XSU(IAC1) - BX*XSU(IAC2)
  149.           end if
  150.           dqs2 = 2.D0*YSU(IAC1)-QSI(2)-YSU(IAC2)
  151.           if (abs(dqs2).GT.1.D-10) then
  152.             BY = (QSI(2) - YSU(IAC1))*2.D0/dqs2
  153.             QSI(2) =(1.D0+BY)*YSU(IAC1) - BY*YSU(IAC2)
  154.           endif
  155.           dqs3 = 2.D0*ZSU(IAC1)-QSI(3)-ZSU(IAC2)
  156.           if (abs(dqs3).GT.1.D-10) then
  157.             BZ = (QSI(3) - ZSU(IAC1))*2.D0/dqs3
  158.             QSI(3) =(1.D0+BZ)*ZSU(IAC1) - BZ*ZSU(IAC2)
  159.           endif
  160.           ITACC = ITACC + ITACN
  161.         ENDIF
  162.  
  163.         i_z = IAC2
  164.         IAC2 = IAC1
  165.         IAC1 = IACC
  166.         IACC = i_z
  167.  
  168.         XSU(IACC) = QSI(1)
  169.         YSU(IACC) = QSI(2)
  170.         ZSU(IACC) = QSI(3)
  171.  
  172. C-FIN---- ITERATIONS ----------------------
  173.       GOTO 610
  174.  
  175. C==========  Cas IDI = 2  ==============================================
  176.  700  CONTINUE
  177. C
  178. C     CAS 2 ELEMENTS SURFACIQUES
  179. C
  180.       KELT = MELE
  181. * Cas particulier des POLYGONES : meleme.itypel = 32 et kelt = 108 + nbnoe
  182.       IF (MELE.EQ.32) THEN
  183.         KELT = 108 + NBMA1
  184.       ENDIF
  185.       EPSIp2 = EPSI*EPSI
  186.  
  187. C-------------------- DEBUT DES ITERATIONS ------------------------
  188.  710  CONTINUE
  189.         ITER = ITER + 1
  190.         IF (ITER.GT.100) THEN
  191.           IROUT = ITER
  192.           GOTO 900
  193.         ENDIF
  194.  
  195. * Cas particulier des POLYGONEs
  196.         IF (MELE.EQ.32) THEN
  197.           CALL SHPOLY(QSI(1),QSI(2),QSI(3),KELT,SHP,iret)
  198.         ELSE
  199.           CALL SHAPE(QSI(1),QSI(2),QSI(3),MELE,SHP,iret)
  200.         ENDIF
  201.         if (iret.eq.0) then
  202.           MOTERR(1:4)=NOMS(MELE)
  203.           CALL ERREUR(68)
  204.           IROUT = 9999
  205.           GOTO 900
  206.         endif
  207.  
  208.         F_1 = XZER
  209.         F_2 = XZER
  210.         DO ik = 1, NBMA1
  211.           F_1 = F_1 + SHP(1,ik) * XEL(1,ik)
  212.           F_2 = F_2 + SHP(1,ik) * XEL(2,ik)
  213.         ENDDO
  214.         F_1 = F_1 - XPG(1)
  215.         F_2 = F_2 - XPG(2)
  216.  
  217. c*        F_DIST = (F_1 * F_1) + (F_2 * F_2)
  218. c*        IF (F_DIST .LT. EPSIp2) GOTO 900
  219.  
  220.         DXDQSI = XZER
  221.         DXDETA = XZER
  222.         DYDQSI = XZER
  223.         DYDETA = XZER
  224.         DO ik = 1, NBMA1
  225.           DXDQSI = DXDQSI + SHP(2,ik) * XEL(1,ik)
  226.           DXDETA = DXDETA + SHP(3,ik) * XEL(1,ik)
  227.           DYDQSI = DYDQSI + SHP(2,ik) * XEL(2,ik)
  228.           DYDETA = DYDETA + SHP(3,ik) * XEL(2,ik)
  229.         ENDDO
  230.         DJAC=DXDQSI*DYDETA-DXDETA*DYDQSI
  231.         XX = DJAC
  232.         IF (DJAC.NE.XZER) XX = XUN/DJAC
  233.         DINV(1,1) =  DYDETA*XX
  234.         DINV(2,1) = -DXDETA*XX
  235.         DINV(1,2) = -DYDQSI*XX
  236.         DINV(2,2) =  DXDQSI*XX
  237.  
  238.         dqs1 = DINV(1,1)*F_1 + DINV(2,1)*F_2
  239.         dqs2 = DINV(1,2)*F_1 + DINV(2,2)*F_2
  240.         QSI(1) = QSI(1) - dqs1
  241.         QSI(2) = QSI(2) - dqs2
  242.  
  243.         dDQ = (dqs1 * dqs1) + (dqs2 * dqs2)
  244.         IF (dDQ.LT.EPSIp2) GOTO 900
  245.  
  246. C Acceleration de convergence :
  247.         IF (ITER.EQ.ITACC) THEN
  248.           dqs1 = 2.D0*XSU(IAC1)-QSI(1)-XSU(IAC2)
  249.           if (abs(dqs1).GT.1.D-10) then
  250.             BX = (QSI(1) - XSU(IAC1))*2.D0/dqs1
  251.             QSI(1) =(1.D0+BX)*XSU(IAC1) - BX*XSU(IAC2)
  252.           end if
  253.           dqs2 = 2.D0*YSU(IAC1)-QSI(2)-YSU(IAC2)
  254.           if (abs(dqs2).GT.1.D-10) then
  255.             BY = (QSI(2) - YSU(IAC1))*2.D0/dqs2
  256.             QSI(2) =(1.D0+BY)*YSU(IAC1) - BY*YSU(IAC2)
  257.           endif
  258.           ITACC = ITACC + ITACN
  259.         ENDIF
  260.  
  261.         i_z = IAC2
  262.         IAC2 = IAC1
  263.         IAC1 = IACC
  264.         IACC = i_z
  265.  
  266.         XSU(IACC) = QSI(1)
  267.         YSU(IACC) = QSI(2)
  268.  
  269. C-FIN---- ITERATIONS ----------------------
  270.       GOTO 710
  271.  
  272. C==========  Cas IDI = 1  ==============================================
  273.  800  CONTINUE
  274. C Methode de Newton-Raphson pour resoudre X(QSI)-XPG(1)=0
  275.  810  CONTINUE
  276.         ITER = ITER+1
  277.         IF (ITER.GT.100) THEN
  278.           IROUT = ITER
  279.           GOTO 900
  280.         ENDIF
  281.  
  282.         CALL SHAPE(QSI(1),QSI(2),QSI(3),MELE,SHP,iret)
  283.         if (iret.eq.0) then
  284.           MOTERR(1:4)=NOMS(MELE)
  285.           CALL ERREUR(68)
  286.           IROUT = 9999
  287.           GOTO 900
  288.         endif
  289.  
  290.         F_1   =XZER
  291.         dXdQSI=XZER
  292.         DO ik = 1, NBMA1
  293.           F_1    = F_1    + SHP(1,ik) * XEL(1,ik)
  294.           dXdQSI = dXdQSI + SHP(2,ik) * XEL(1,ik)
  295.         ENDDO
  296.         F_1 = F_1 - XPG(1)
  297. c*        F_DIST = ABS(F_1)
  298. c*        IF (F_DIST .LT. EPSIp2) GOTO 900
  299.         IF (dXdQSI.NE.XZER) THEN
  300.           dQS1 = F_1 / dXdQSI
  301.           QSI(1)=QSI(1)-dQS1
  302.         ELSE
  303.           dQS1=XZER
  304.         ENDIF
  305. C Test de convergence
  306.         IF (ABS(dQS1).LT.EPSI) GOTO 900
  307. C Acceleration de convergence :
  308.         IF (ITER.EQ.ITACC) THEN
  309.           dqs1 = 2.D0*XSU(IAC1)-QSI(1)-XSU(IAC2)
  310.           if (abs(dqs1).GT.1.D-10) then
  311.             BX = (QSI(1) - XSU(IAC1))*2.D0/dqs1
  312.             QSI(1) =(1.D0+BX)*XSU(IAC1) - BX*XSU(IAC2)
  313.           end if
  314.           ITACC = ITACC + ITACN
  315.         ENDIF
  316.  
  317.         i_z = IAC2
  318.         IAC2 = IAC1
  319.         IAC1 = IACC
  320.         IACC = i_z
  321.  
  322.         XSU(IACC) = QSI(1)
  323.  
  324. C-FIN---- ITERATIONS ----------------------
  325.       GOTO 810
  326.  
  327. C========  Sortie  =============
  328.  900  CONTINUE
  329. c*      write(6,*) 'QSIJS',IROUT,ITER,QSI(1),QSI(2),QSI(3),F_DIST
  330.  
  331. c      return
  332.       END
  333.  
  334.  
  335.  

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