Télécharger la plaquette Cast3M


fr en

Télécharger act3.eso

Retour à la liste

Numérotation des lignes : 

act3

  1. C ACT3      SOURCE    PV        22/02/24    21:15:01     11297          
  2.       SUBROUTINE ACT3
  3. C--------------------------------------------------------------------
  4. C              ACCELERATION POUR DES CHAMPOINTS
  5. C
  6. C   APPELE DANS INCREME
  7. C
  8. C   LIT DANS L'ORDRE 3 CHAMPS PRIMAUX
  9. C                 ET 4 CHAMPS DE DUAUX
  10. C
  11. C   REND UN CHAMP PRIMAL EXTRAPOLE QUI MINIMISE LE CHAMP DUAL.
  12. C   METHODE : ON ESTIME L'APPLICATION TANGENTE ET ON L'UTILISE
  13. C   POUR RESOUDRE LE PROBLEME
  14. C
  15. C--------------------------------------------------------------------
  16. C
  17. C   AOUT 97 : RAJOUTE L'ACCELERATION AVEC 2 OU 3 CHAMPS DE DEPLACEMENTS
  18. C
  19. C   CHAT 05/10/2007 : suppression de la limitation à certaines composantes
  20. C
  21. C   PM 19/07/2007 : cas à 2 points mettait IDCH0 à 0 au lieu de IDCH1
  22. C
  23. C--------------------------------------------------------------------
  24.       IMPLICIT INTEGER(I-N)
  25.       IMPLICIT REAL*8(A-H,O-Z)
  26. -INC CCREEL
  27. -INC SMCHPOI
  28. -INC SMELEME
  29. -INC SMLMOTS
  30.  
  31. -INC PPARAM
  32. -INC CCOPTIO
  33.         xref=xszpre
  34. C
  35. C  FABRICATION DES LISTMOTS
  36. C
  37.       JGN=LOCOMP
  38.       JGM=100
  39.       SEGINI MLMOTS
  40.       JGM=0
  41.  
  42.       MLMOTF=MLMOTS
  43. *
  44.       CALL LIROBJ('CHPOINT ',ICHD2,1,IRETOU)
  45.       CALL LIROBJ('CHPOINT ',ICHD1,1,IRETOU)
  46.       CALL LIROBJ('CHPOINT ',ICHD0,1,IRETOU)
  47.       CALL LIROBJ('CHPOINT ',IFO3 ,1,IRETOU)
  48.       CALL LIROBJ('CHPOINT ',IFO2 ,1,IRETOU)
  49.       CALL LIROBJ('CHPOINT ',IFO1 ,1,IRETOU)
  50.       CALL LIROBJ('CHPOINT ',IFO0 ,1,IRETOU)
  51.  
  52.       CALL ACTOBJ('CHPOINT ',ICHD2,1)
  53.       CALL ACTOBJ('CHPOINT ',ICHD1,1)
  54.       CALL ACTOBJ('CHPOINT ',ICHD0,1)
  55.       CALL ACTOBJ('CHPOINT ',IFO3 ,1)
  56.       CALL ACTOBJ('CHPOINT ',IFO2 ,1)
  57.       CALL ACTOBJ('CHPOINT ',IFO1 ,1)
  58.       CALL ACTOBJ('CHPOINT ',IFO0 ,1)
  59.       IF (IERR.NE.0) RETURN
  60.  
  61. *     Acceleration avec 3 points
  62.       IF (ICHD2.EQ.ICHD1) THEN
  63.        write(6,*) ' acceleration sur 3 points '
  64.           ICHD2 = 0
  65.           IFO3  = 0
  66.        ENDIF
  67.  
  68. *     Acceleration avec 2 points
  69.       IF (ICHD1.EQ.ICHD0) THEN
  70.        write(6,*) ' acceleration sur 2 points '
  71.           ICHD2 = 0
  72. *PM          ICHD0 = 0
  73.           ICHD1 = 0
  74.           IFO3  = 0
  75.           IFO2  = 0
  76.        ENDIF
  77.  
  78. *
  79. *     Recherche des noms de composante des forces
  80. *
  81.       MCHPOI=IFO3
  82.       IPA=1
  83.     1 CONTINUE
  84.       IF (MCHPOI.NE.0) THEN
  85.          SEGACT MCHPOI
  86.          DO IA=1,IPCHP(/1)
  87.             MSOUPO=IPCHP(IA)
  88.             SEGACT MSOUPO
  89.             DO IB=1,NOCOMP(/2)
  90.                DO IC=1,JGM
  91.                   IF (MOTS(IC).EQ.NOCOMP(IB)) GO TO 2
  92.                ENDDO
  93. **             write (6,*) ' act3 ncomp ',nocomp(ib)
  94.                if (nocomp(ib).eq.'FLX ') goto 2
  95.  
  96.  
  97. *         on ajoute la nouvelle composante à la liste
  98.                JGM=JGM+1
  99.                IF (JGM.GT.MOTS(/2)) THEN
  100.                   SEGADJ MLMOTS
  101.                ENDIF
  102.                MOTS(JGM)=NOCOMP(IB)
  103.  2             CONTINUE
  104.             ENDDO
  105.          ENDDO
  106.       ENDIF
  107. *     On fait de même pour les trois autres champs en entrée
  108.       IF (IPA.EQ.1) THEN
  109.           IPA=2
  110.           MCHPOI=IFO2
  111.           GO TO 1
  112.       ELSEIF(IPA.EQ.2) THEN
  113.           IPA=3
  114.           MCHPOI=IFO1
  115.           GO TO 1
  116.       ELSEIF(IPA.EQ.3) THEN
  117.           IPA=4
  118.           MCHPOI=IFO0
  119.           GO TO 1
  120.       ENDIF
  121.       SEGADJ MLMOTS
  122. *
  123. *     On utilise un procédé d'orthogonalisation de Schmidt
  124. *     pour obtenir une base de l'espace engendrée par les quatre forces
  125. *     puis on effectue la projection (de coordonnées XL, XM et XN) du
  126. *     champ de force nulle sur cet espace
  127. *
  128.       IF (IFO3.NE.0) THEN
  129.  
  130. *     CAS  A 4 POINTS :
  131. *     -----------------
  132. *     famille des 3 vecteurs variations du champ dual
  133.       CALL ADCHPO(IFO1,IFO0,IFO10,1D0,-1D0)
  134.       CALL ADCHPO(IFO2,IFO0,IFO20,1D0,-1D0)
  135.       CALL ADCHPO(IFO3,IFO0,IFO30,1D0,-1D0)
  136. *
  137.       CALL XTX1(IFO10,F12)
  138.       CALL XTX1(IFO20,F22)
  139.       CALL XTX1(IFO30,F32)
  140.       if(f12.LT.XPETIT*1d2) f12=XPETIT*1d2
  141.       if(f22.LT.XPETIT*1d2) f22=XPETIT*1d2
  142.       if(f32.LT.XPETIT*1d2) f32=XPETIT*1d2
  143. *     write (6,*) ' dans act3 f12/f22 ',f12/f22,f12/f32
  144. *
  145. *     Détection de flip-flop ?
  146. *     (1 et 3 sont trop proches, malgré 2 différent de 1)
  147.       if (f12/f22.gt.1.2d0 .and. f12/f32.gt.0.707d0 .and.
  148.      >       f12/f32.lt.1.414d0) then
  149. *       act3 : plan B
  150.         CALL ERREUR(858)
  151.         XL = (sqrt(f12)+sqrt(f22))/(2*sqrt(f12))
  152.         XM = 0.
  153.         XN = 0.
  154.         GOTO 1010
  155.       endif
  156. *
  157. *     Construction d'une famille libre orthogonale
  158.       CALL XTY1(IFO10,IFO20,MLMOTF,MLMOTF,F1F2)
  159.       CALL XTY1(IFO10,IFO30,MLMOTF,MLMOTF,F1F3)
  160.       CALL XTY1(IFO20,IFO30,MLMOTF,MLMOTF,F2F3)
  161.       CALL XTY1(IFO0 ,IFO10,MLMOTF,MLMOTF,F0F1)
  162.       CALL XTY1(IFO0 ,IFO20,MLMOTF,MLMOTF,F0F2)
  163.       CALL XTY1(IFO0 ,IFO30,MLMOTF,MLMOTF,F0F3)
  164. *
  165. *  normalisation pour ameliorer la precision
  166. *  les iteres successives pouvant diminer pas mal si on converge
  167. *     A1  =  F12
  168. *     B1  =  F1F2
  169. *     C1  =  F1F3
  170. *     XJ1 = -F0F1
  171. *     D1  =  F1F2
  172. *     E1  =  F22
  173. *     F1  =  F2F3
  174. *     XK1 = -F0F2
  175. *     G1  =  F1F3
  176. *     H1  =  F2F3
  177. *     XI1 =  F32
  178. *     XL1 = -F0F3
  179.       A1  =  1.d0
  180.       B1  =  F1F2/f12
  181.       C1  =  F1F3/f12
  182.       XJ1 = -F0F1/f12
  183.       D1  =  F1F2/f22
  184.       E1  =  1.D0
  185.       F1  =  F2F3/f22
  186.       XK1 = -F0F2/f22
  187.       G1  =  F1F3/f32
  188.       H1  =  F2F3/f32
  189.       XI1 =  1.d0
  190.       XL1 = -F0F3/f32
  191.  
  192.  
  193.  
  194. *     write (6,*) ' f12 f22 f32 ',f12,f22,f32
  195.  
  196. *     Calcul des coordonnées de IFO0 dans la famille libre
  197. *     et élimination des vecteurs liés
  198.       IF (ABS(XI1) .LT. xpetit) GOTO 2000
  199. *
  200.       A2  = A1  - ( G1*C1/XI1)
  201.       B2  = B1  - ( H1*C1/XI1)
  202.       XJ2 = XJ1 - (XL1*C1/XI1)
  203.       D2  = D1  - ( G1*F1/XI1)
  204.       E2  = E1  - ( H1*F1/XI1)
  205.       XK2 = XK1 - (XL1*F1/XI1)
  206. *
  207.       IF (E2 .LE. e1*xref) THEN
  208. *           write (6,*) ' act3 e2 e1 ',e2,e1
  209.          GOTO 2000
  210.       ENDIF
  211. *
  212.       A3  = A2 -( D2*B2/E2)
  213.       XJ3 = XJ2-(XK2*B2/E2)
  214. *
  215.       IF (A3 .LE. a1*xref) THEN
  216. *           write (6,*) ' act3 a3 a1 ',a3,a1
  217.          GOTO 2000
  218.       ENDIF
  219.       IF (ABS(A2) .LE. A1*1d-6) THEN
  220. *           write (6,*) ' act3 a2 a1 ',a2,a1
  221.          GOTO 2000
  222.       ENDIF
  223.  
  224.       XL=XJ3/A3
  225.       if (xl .gt.1.D0) then
  226. **         write(6,*) ' act3 xl > 0 ==> pas d acceleration'
  227. **         goto 2000
  228.       endif
  229. *
  230.       B4  = B2
  231.       XJ4 = XJ2 - (A2*XL)
  232.  
  233.  
  234. *
  235.       IF (ABS(B4) .LE. ABS(B1)*1D-6) THEN
  236. *           write (6,*) ' act3 b4 b1 ',b4,b1
  237.          GOTO 2000
  238.       ENDIF
  239. *
  240.       XM  = XJ4/B4
  241. *
  242.       C5  = C1
  243.       XJ5 = XJ1 - (A1*XL) - (B1*XM)
  244. *
  245.       IF (ABS(C5) .LT. xref) GOTO 2000
  246. *
  247.       XN=XJ5/C5
  248. *
  249.       ref=sqrt(a1**2+b1**2+c1**2)
  250. *     write (6,*) ' act3 ',ref,XI1,A3,b4,c5
  251. *     write (6,*) ' act3 a3 e2 xi1 ',a3,e2,xi1
  252.       if (a3.lt.a1*xref.or.e2.lt.e1*xref) then
  253.       write (6,*) ' act3 operateur secant non positif'
  254.        goto 2000
  255.       endif
  256.       GOTO 1010
  257.  
  258.  2000 CONTINUE
  259. *     le systeme est singulier on essaye avec un vecteur de moins
  260. *
  261. *     write (6,*) ' act3 ',xi1,e2,e1,a3,a1,b4,b1,c5
  262.       IF (ABS(E1) .LT. xref) GOTO 2100
  263. *
  264.       A2  = A1  - (D1 *B1/E1)
  265.       XJ2 = XJ1 - (XK1*B1/E1)
  266. *
  267.       IF (A2 .LE. A1*xref) THEN
  268. *        write (6,*) ' act3 a2 a1 ',a2,a1
  269.          GOTO 2100
  270.       ENDIF
  271. *
  272.       XL  = XJ2/A2
  273.       if (xl .gt.1.D0) then
  274. **         write(6,*) ' act3 xl > 0 ==> pas d acceleration'
  275. **         goto 2100
  276.       endif
  277. *
  278.       B3  = B1
  279.       XJ3 = XJ1 - (A1*XL)
  280. *
  281.       IF (ABS(B3) .LT. xref) GOTO 2100
  282. *
  283.       XM = XJ3/B3
  284. *
  285.       XN = 0
  286.  
  287.       if (a2.lt.a1*xref) then
  288. *     write (6,*) ' act3 a2 a1  ',a2,e1
  289.        write (6,*) ' act3 operateur secant non positif'
  290.        goto 2100
  291.       endif
  292. *     act3 : reduction a %i1 dimensions
  293.       INTERR(1)=2
  294.        CALL ERREUR(859)
  295.       GOTO 1010
  296.  
  297.  2100 CONTINUE
  298. *     le systeme est singulier on essaye avec un vecteur de moins
  299. *
  300.       IF (ABS(A1) .LT. xref ) GOTO 1000
  301. *
  302.       XL  = XJ1/A1
  303.       if (xl .gt.1.D0) then
  304. **         write(6,*) ' act3 xl > 0 ==> pas d acceleration'
  305. **         goto 1000
  306.       endif
  307. *
  308.       XM  = 0
  309.       XN  = 0
  310.  
  311. *     act3 : reduction a %i1 dimensions
  312. *     write (6,*) ' act3 a1 ',a1
  313.       INTERR(1)=1
  314.  
  315.       CALL ERREUR(859)
  316.       GOTO 1010
  317.  
  318.  1000 CONTINUE
  319. *     act3 : pas d accélération
  320.       call erreur(860)
  321.       XL = 0
  322.       XM = 0
  323.       XN = 0
  324.  
  325.  1010 CONTINUE
  326. *
  327. *     WRITE (6,*) ' f12*f22/(f1f2*2) ',f12*f22/(f1f2**2)
  328. *     WRITE (6,*) ' f12*f32/(f1f3*2) ',f12*f32/(f1f3**2)
  329. *     WRITE (6,*) ' f22*f32/(f2f3*2) ',f22*f32/(f2f3**2)
  330. *     write (6,*) ' xi1,e2,a3,b4,c5 ',xi1,e2,a3,b4,c5
  331. *     write (6,*) ' coefficients acc ',-xl,-xm,-xn
  332. *     xl=sign(min(abs(xl),1d1),xl)
  333. *     xm=sign(min(abs(xm),1d1),xm)
  334. *     xn=sign(min(abs(xn),1d1),xn)
  335. *     write (6,*) ' coefficients cor ',-xl,-xm,-xn
  336.  
  337. *     Vérification qu'on ne retombe pas sur l'un des duals entrés
  338. *     Sinon, on réessaye en tenant compte de moins de vecteurs
  339.       if (abs(xl-1).lt.1e-1 .and. abs(xm).lt.1e-1
  340.      >                      .and. abs(xn).lt.1e-1) goto 1000
  341.       if (abs(xm-1).lt.1e-1 .and. abs(xn).lt.1e-1
  342.      >                      .and. abs(xl).lt.1e-1) goto 2100
  343.       if (abs(xn-1).lt.1e-1 .and. abs(xl).lt.1e-1
  344.      >                      .and. abs(xm).lt.1e-1) goto 2000
  345. *
  346.       CALL ADCHPO(ICHD0,ICHD1,ITEMP4,XL+XM+XN,XM+XN)
  347.       CALL ADCHPO(ICHD2,ITEMP4,ITEMP5,XN,1D0)
  348. *
  349. *     ajustement de la nature du champ
  350. *     on prend la meme nature que ICHDO
  351.       MCHPOI = ICHD0
  352.       MCHPO1 = ITEMP5
  353.       SEGACT MCHPOI,MCHPO1
  354.       NAT    = JATTRI(/1)
  355.       NSOUPO = MCHPO1.IPCHP(/1)
  356.       SEGADJ MCHPO1
  357.       DO 1020 I=1,NAT
  358.         MCHPO1.JATTRI(I)=JATTRI(I)
  359.  1020 CONTINUE
  360. *
  361. *     Affectation
  362. *
  363.       CALL ACTOBJ('CHPOINT ',ITEMP5,1)
  364.       CALL ECROBJ('CHPOINT ',ITEMP5)
  365. *
  366. *     Quelques destructions
  367. *
  368.       MLMOTS=MLMOTF
  369.       SEGSUP MLMOTS
  370.       CALL DTCHPO(IFO10)
  371.       CALL DTCHPO(IFO20)
  372.       CALL DTCHPO(IFO30)
  373.       CALL DTCHPO(ITEMP4)
  374.       RETURN
  375.       ENDIF
  376.  
  377. *     CAS  A 3 POINTS :
  378. *     -----------------
  379. *     famille des 2 vecteurs variations du champ dual
  380.       IF (IFO2.NE.0.AND.IFO3.EQ.0) THEN
  381.       CALL ADCHPO(IFO1,IFO0,IFO10,1D0,-1D0)
  382.       CALL ADCHPO(IFO2,IFO0,IFO20,1D0,-1D0)
  383. *
  384.       CALL XTX1(IFO10,F12)
  385.       CALL XTX1(IFO20,F22)
  386. *
  387. *     Construction d'une famille libre orthogonale
  388.       CALL XTY1(IFO10,IFO20,MLMOTF,MLMOTF,F1F2)
  389.       CALL XTY1(IFO0 ,IFO10,MLMOTF,MLMOTF,F0F1)
  390.       CALL XTY1(IFO0 ,IFO20,MLMOTF,MLMOTF,F0F2)
  391. *
  392. *      A1  =  F12
  393. *      B1  =  F1F2
  394. *      XJ1 = -F0F1
  395. *      D1  =  F1F2
  396. *      E1  =  F22
  397. *      XK1 = -F0F2
  398.       A1  =  1.d0
  399.       B1  =  F1F2/F12
  400.       XJ1 = -F0F1/F12
  401.       D1  =  F1F2/F22
  402.       E1  =  1.D0
  403.       XK1 = -F0F2/F22
  404. *
  405. *     Calcul des coordonnées de IFO0 dans la famille libre
  406. *     et élimination des vecteurs liés
  407.       IF (ABS(E1) .LT. xref ) GOTO 2200
  408. *
  409.       A2  = A1  - (D1 *B1/E1)
  410.       XJ2 = XJ1 - (XK1*B1/E1)
  411. *
  412.       IF (A2 .LE. A1*xref) THEN
  413. *        write (6,*) ' act3 a2 a1 ',a2,a1
  414.          GOTO 2200
  415.       ENDIF
  416. *
  417.       XL  = XJ2/A2
  418. *
  419.       B3  = B1
  420.       XJ3 = XJ1 - (A1*XL)
  421. *
  422.       IF (ABS(B3) .LT. xref) GOTO 2200
  423. *
  424.       XM  = XJ3/B3
  425. *
  426.       XN  = 0
  427. *
  428. *     write (6,*) ' act3 xl xm ',xl,xm
  429.       GOTO 1210
  430.  
  431.  2200 CONTINUE
  432. *     le systeme est singulier on essaye avec un vecteur de moins
  433. *
  434.       IF (A1 .LT. xref ) GOTO 2300
  435. *
  436.       XL  = XJ1/A1
  437. *
  438.       XM  = 0
  439.       XN  = 0
  440. *
  441. *     act3 : reduction a %i1 dimensions
  442.       INTERR(1)=1
  443. *      CALL ERREUR(859)
  444.       GOTO 1210
  445.  
  446.  2300 CONTINUE
  447. *     act3 : pas d accélération
  448.       call erreur(860)
  449.       XL  = 0
  450.       XM  = 0
  451.       XN  = 0
  452.  
  453.  1210 CONTINUE
  454. *     Vérification qu'on ne retombe pas sur l'un des duals entrés
  455. *     Sinon, on réessaye en tenant compte de moins de vecteurs
  456.       if (abs(xl-1).lt.1e-1.and.abs(xm).lt.1e-1.and.
  457.      >    abs(xn).lt.1e-1) goto 2300
  458.       if (abs(xm-1).lt.1e-1.and.abs(xn).lt.1e-1.and.
  459.      >    abs(xl).lt.1e-1) goto 2200
  460. *      if (abs(xn-1).lt.1e-1.and.abs(xl).lt.1e-1.and.
  461. *     >    abs(xm).lt.1e-1) goto 2300
  462. *
  463.       CALL ADCHPO(ICHD0,ICHD1,ITEMP5,XL+XM,XM)
  464. *
  465. *     ajustement de la nature du champ
  466. *     on prend la meme nature que ICHDO
  467. *
  468.       MCHPOI = ICHD0
  469.       MCHPO1 = ITEMP5
  470.       SEGACT MCHPOI,MCHPO1
  471.       NAT    = JATTRI(/1)
  472.       NSOUPO = MCHPO1.IPCHP(/1)
  473.       SEGADJ MCHPO1
  474.       DO 1220 I=1,NAT
  475.         MCHPO1.JATTRI(I)=JATTRI(I)
  476.  1220 CONTINUE
  477. *
  478. *     Affectation
  479. *
  480.       CALL ACTOBJ('CHPOINT ',ITEMP5,1)
  481.       CALL ECROBJ('CHPOINT ',ITEMP5)
  482. *
  483. *     Quelques destructions
  484. *
  485.       MLMOTS=MLMOTF
  486.       SEGSUP MLMOTS
  487.       CALL DTCHPO(IFO10)
  488.       CALL DTCHPO(IFO20)
  489.       RETURN
  490.       ENDIF
  491.  
  492. *     CAS  A 2 POINTS :
  493. *     -----------------
  494.       IF (IFO2.EQ.0) THEN
  495.       CALL ADCHPO(IFO1,IFO0,IFO10,1D0,-1D0)
  496. *
  497.       CALL XTX1(IFO10,F12)
  498. *
  499. *     Un seul vecteur = famille libre d'emblée
  500.       CALL XTY1(IFO0,IFO10,MLMOTF,MLMOTF,F0F1)
  501. *
  502. *      A1  =  F12
  503. *      XJ1 = -F0F1
  504.       A1  =  1.D0
  505.       XJ1 = -F0F1/F12
  506. *
  507.       IF (ABS(A1) .LT. xref ) GOTO 1100
  508. *
  509.       XL  = XJ1/A1
  510. *
  511.       XM  = 0
  512.       XN  = 0
  513.       GOTO 1110
  514.  
  515.  1100 CONTINUE
  516. *     act3 : pas d accélération
  517.       call erreur(860)
  518.       XL = 0
  519.       XM = 0
  520.       XN = 0
  521.  
  522.  1110 CONTINUE
  523.       if (abs(xl-1).lt.1e-1) goto 1100
  524.       CALL MUCHPO(ICHD0,XL,ITEMP5,1)
  525. *     write (6,*) ' act3 xl vaut ',xl
  526. *
  527. *     ajustement de la nature du champ
  528. *     on prend la meme nature que ICHDO
  529. *
  530.       MCHPOI = ICHD0
  531.       MCHPO1 = ITEMP5
  532.       SEGACT MCHPOI,MCHPO1
  533.       NAT = JATTRI(/1)
  534.       NSOUPO = MCHPO1.IPCHP(/1)
  535.       SEGADJ MCHPO1
  536.       DO 1120 I=1,NAT
  537.         MCHPO1.JATTRI(I)=JATTRI(I)
  538.  1120 CONTINUE
  539. *
  540. *     Affectation
  541. *
  542.       CALL ACTOBJ('CHPOINT ',ITEMP5,1)
  543.       CALL ECROBJ('CHPOINT ',ITEMP5)
  544. *
  545. *     QUELQUES DESTRUCTIONS
  546. *
  547.       MLMOTS=MLMOTF
  548.       SEGSUP MLMOTS
  549.       CALL DTCHPO(IFO10)
  550.       ENDIF
  551.       END
  552.  
  553.  
  554.  
  555.  
  556.  
  557.  
  558.  

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