Télécharger la plaquette Cast3M


fr en

Télécharger mondes.eso

Retour à la liste

Numérotation des lignes : 

mondes

  1. C MONDES    SOURCE    PV090527  26/02/20    21:15:01     12479          
  2.       SUBROUTINE MONDES(MMATRX,MVECTX,NOEN,ISOUCI,lagdua)
  3. C
  4. C  **** EXECUTE LA SOLUTION X DE  (Lt D L) X=F
  5. C
  6. CMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMB
  7. CMB
  8. CMB   Plutot la solution de L.D.Lt ou L.D.Mt (cas non symétrique)
  9. CMB   Elle devrait dons s'appeller DESMON et non MONDES.
  10. CMB
  11. CMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMB
  12. C
  13.       IMPLICIT INTEGER(I-N)
  14.       IMPLICIT REAL*8(A-H,O-Z)
  15.       PARAMETER (LPCRAY=10000)
  16.       INTEGER OOOVAL
  17.  
  18. C ... Variable décrivant l'EXIStence du Triangle Supérieur différent
  19. C     de l'inférieur (cas non symétrique) ...
  20.       LOGICAL EXISTS
  21.  
  22. -INC SMMATRI
  23. -INC SMELEME
  24. -INC SMVECTD
  25.       POINTEUR MVECT4.MVECTD
  26.  
  27. -INC PPARAM
  28. -INC CCOPTIO
  29. -INC CCREEL
  30. -INC SILICRE
  31. -INC CCASSIS
  32. -INC TMTRAV
  33.  
  34.       SEGMENT ITTRV(MULRE)
  35.       SEGMENT,ITRAA(NENSL)
  36.       SEGMENT/BID/(IBIDON(IIMAX))
  37.       logical toutmem
  38.       mvect4 = 0
  39.       mtrav=0
  40.       toutmem=.true.
  41.       DNORMA=0.D0
  42.       DNORMB=0.D0
  43.        crit=xpetit
  44.        crit2 = 0.D0
  45.       call oooprl(1)
  46.       nbthr=nbthrs
  47.       CALL OOOMRU(1)
  48.       IVALMA=0
  49.       IF(IIMPI.EQ.3) WRITE(IOIMP,1000)MMATRX,MVECTX
  50.  1000 FORMAT('  SUBROUTINE  MONDES : POINTEUR DE LA MATRICE=',I5,
  51.      1                          '    POINTEUR DU VECTEUR=',I5)
  52. C
  53. C  ****  ACTIVATION DES SEGMENTS
  54. C
  55.       MMATRI=MMATRX
  56.       SEGACT,MMATRI*MOD
  57. *  BEC00SE OPTIMISEUR
  58.       ITRAA=MMATRI
  59.       NENSL=0
  60.       IF(NENS.NE.0) THEN
  61.        NENSL=NENS
  62.        SEGINI,ITRAA
  63.       ENDIF
  64.       NE1=NENSL
  65.  
  66.       EXISTS=.FALSE.
  67.       IF(IILIGS.NE.0) THEN
  68.          MILIGN=IILIGS
  69.          SEGACT MILIGN
  70.          IF(ILIGN(/1).GT.0) EXISTS=.TRUE.
  71.       ENDIF
  72.  
  73.       MILIGN=IILIGN
  74.       SEGACT,MILIGN
  75.       INC=IPNO(/1)
  76.       nbthr=min(nbthrs,inc/1000+1)
  77.  
  78.       MVECTD=MVECTX
  79.       SEGACT MVECTD*MOD
  80.  
  81.       MDNOR=IDNORM
  82.       SEGACT MDNOR
  83.       IF(IDNORD.GT.0) THEN
  84.          MDNO1=IDNORD
  85.          SEGACT MDNO1
  86.       ELSE
  87.          MDNO1=MDNOR
  88.       ENDIF
  89.  
  90.       MDIAG=IDIAG
  91.  
  92. C ... MULRE = nombre de seconds membres ...
  93.       MULRE = VECTBB(/1)/INC
  94. C ... MUINC ne servira que comme une borne sur les indices de boucles ...
  95.       MUINC = ( MULRE-1)*INC
  96.       SEGINI ITTRV
  97.       inc=vectbb(/1)
  98. *  mvect1 sauve les forces
  99. *  mvect2 sauve le deplacement
  100. *  mvect3 vecteur de travail
  101. *  mvectd evolue: force a resoudre (residu) puis increment de deplacement
  102.       if (mulre.eq.1.and..not.exists) segini mvect1,mvect2
  103.       INC=IPNO(/1)
  104.  
  105. C ... Multiplication du second membre par les DNOR(.) (à cause du
  106. C     conditionnement des matrices de blocages) ...
  107.       DNORMA=0.D0
  108.       DO 450 K=0,MUINC,INC
  109.          DO 45 I=1,INC
  110.             VECTBB(I+K)=VECTBB(I+K)*MDNO1.DNOR(I)
  111. **          DNORMA=MAX(DNORMA,ABS(VECTBB(I+K)))
  112.  45      CONTINUE
  113.  450  CONTINUE
  114. **    DNORMA=DNORMA*xzprec*xzprec
  115.       DNORMA=xpetit/xzprec
  116.  
  117. C ... Détection du premier terme important du second membre ...
  118.       II=0
  119.       DO 451 K=0,MUINC,INC
  120.          II=II+1
  121.          DO 452 I=1,INC
  122.             IF(ABS(VECTBB(I+K)).LT.DNORMA) GO TO 452
  123. C       ... Le numéro du noeud concerné par le premier terme important va à ITTRV ...
  124.             ITTRV(II)=IPNO(I)
  125.             GO TO 451
  126.   452    CONTINUE
  127.   451 CONTINUE
  128.  
  129. C ... NNOE = nombre de noeuds concernés ...
  130.       NNOE=ILIGN(/1)
  131.       IDEP=NNOE
  132. C ... On cherche le minimum (non nul) des ITTRV que l'on met dans IDEP ...
  133.       DO 453 I=1,MULRE
  134.          IF(ITTRV(I).LT.IDEP) THEN
  135.             IF(ITTRV(I).NE.0) THEN
  136.                IDEP=ITTRV(I)
  137.             ELSE
  138.                IDEP=1
  139.                GO TO 4530
  140.             ENDIF
  141.          ENDIF
  142.  453  CONTINUE
  143.  4530 CONTINUE
  144.       SEGDES MDNOR
  145. *
  146. *   debut de la boucle d'amelioration de la precision si necessaire
  147. *
  148. *   sauver le champs de forces initial
  149.       if(mulre.eq.1.and..not.exists) then
  150.        do i=1,inc
  151.          mvect1.vectbb(i)=vectbb(i)
  152.        enddo
  153.        icorrec=0
  154.       endif
  155. 10000 continue
  156.       IAA=0
  157.  
  158.  
  159.  
  160.  
  161.  
  162.       SEGACT,MDIAG
  163. C
  164. C  ****  DESCENTE: ON RESOU L*C=B.  EN FAIT ON STOCKE C DANS B
  165. C
  166.       LTRK=MAX(1+LPCRAY,OOOVAL(1,4))
  167.       IIMAX=(((IJMAX +LTRK)/LTRK)+1)*LTRK
  168.       iimax=oooval(1,1)/10
  169. *  test si la place disponible permet de tout mettre en memoire
  170.       xplds=oooval(1,1)-oooval(2,3)
  171.       if (real(nnoe)*ijmax.lt.xplds) iimax=0
  172.       SEGINI /err=3500/ BID
  173.       goto 3501
  174.  3500 continue
  175. *  place insuffisante en memoire centrale. On en fait
  176. *  desactivation de la matrice en partant du debut
  177.       do ivs=1,nnoe
  178.        lign=ilign(ivs)
  179.        segdes lign
  180.       enddo
  181.       segini bid
  182.  3501 continue
  183. C ... IDEP : on commence par ce noeud car le second membre est nul pour
  184. C     tous les precedents ...
  185. *  activation de tous les derniers noeuds car ils sont deja en memoire normalement
  186.       do 2001 ivs=nnoe,1,-1
  187.         lign=ilign(ivs)
  188.         segact/err=2010/ lign
  189.  2001 continue
  190.  2010 continue
  191.       ilmin=ivs+1
  192.       DO 610 IVS=IDEP,NNOE
  193.          LIGN=ILIGN(IVS)
  194.          SEGACT /ERR=50/ LIGN
  195.          IVALMA=IVALMA+VAL(/1)
  196. *pvpv         IF (IVALMA.GT.NGMAXY) GOTO 50
  197.          NA=IMMM(/1)
  198.          IPRELL=IPREL
  199.          DO 611 IA=1,NA
  200. C       ... Si l'inconnue présente un mode d'ensemble ...
  201.             IF(IMMM(IA).NE.0) THEN
  202. C          ... On incrémente le compteur et
  203.                IAA=IAA+1
  204. C          ... on met le N° de la ligne dans ITRAA ...
  205.                ITRAA(IAA)=IPRELL+IA-1
  206.             ENDIF
  207.   611    CONTINUE
  208.          segact lign
  209.          CALL MONDE2(ITTRV(1),IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  210.      >      NA,IPREL,MULRE,INC,IVS,LLOM,dnorma)
  211.  610  CONTINUE
  212. C ... Si on n'a pas eu de pb de mémoire on passe par là ...
  213.       SEGSUP BID
  214. C ... ILMAX vaut le dernier noeud qui a pu être traité ...
  215.       ILMAX=NNOE
  216. C ... On va à la division par la diagonale ...
  217.       GOTO 55
  218.  
  219.   50  CONTINUE
  220. C ... Si on est là, c'est à cause des pb de mémoire ...
  221.       SEGSUP BID
  222.       toutmem=.false.
  223. C ... IVS est le N° du LIGN qui n'a pas pu être traité ...
  224. C ... ILMAX = N° du dernier traité ...
  225.       ILMAX=IVS-1
  226. C ... IILMAX = N° du premier à traiter ...
  227.       IILMAX=IVS
  228.       DO 54 IVS=IILMAX,NNOE
  229.   58     CONTINUE
  230.          LIGN=ILIGN(IVS)
  231. C    ... Même si tout à l'heure SEGACT n'a pas marché, maintenant on a supprimé BID ...
  232.          SEGACT /ERR=56/ LIGN
  233.          GOTO 57
  234.   56     CONTINUE
  235. *  EN CAS DE PROBLEME ON FAIT UN PEU DE PLACE
  236.          toutmem=.false.
  237. C    ... Si on a toujours des pb de mémoire, on désactive le LIGN N° ILMAX,
  238.          LIGN=ILIGN(ILMAX)
  239.          if(ilmax.lt.ilmin) SEGDES LIGN*(NOMOD,MRU)
  240. C    ... puis ce dernier est décrémenté ...
  241.          ILMAX=ILMAX-1
  242.          IF (ILMAX.EQ.1) CALL ERREUR(5)
  243. C    ... et on recommence ...
  244.          GOTO 58
  245. C    ... On est là si tout s'est bien passé avec SEGACT ...
  246.   57     CONTINUE
  247.          NA=IMMM(/1)
  248.          IPRELL=IPREL
  249.          DO 612 IA=1,NA
  250.             IF(IMMM(IA).EQ.0) GOTO 612
  251.             IAA=IAA+1
  252.             ITRAA(IAA)=IPRELL+IA-1
  253.   612    CONTINUE
  254.          CALL MONDE2(ITTRV(1),IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  255.      >      NA,IPREL,MULRE,INC,IVS,LLOM,dnorma)
  256.          if (ivs.lt.ilmin) SEGDES,LIGN*(NOMOD,MRU)
  257.    54 CONTINUE
  258.    55 CONTINUE
  259. C
  260. C ... À cet endroit la descente est finie. L'état des LIGN est le suivant :
  261. C ... Ceux de IDEP à ILMAX sont actifs, les autres (si ILMAX < NNOE) sont désactivés ...
  262. C
  263. C  ****  DIVISION PAR LE TERME DIAGONAL  ****
  264. C
  265.       dnorma=0.d0
  266.       idenorm=0
  267.       DO 120 K=0,MUINC,INC
  268.          DO 12 I=1,INC
  269.             J=I+K
  270.             if (2*abs(diag(i)).gt.abs(diag(i))) goto 122
  271.             idenorm=1
  272.  122        continue
  273.             VECTBB(J)=VECTBB(J)*DIAG(I)
  274. **          dnorma=max(dnorma,abs(vectbb(j)))
  275.    12    CONTINUE
  276.  
  277.   120 CONTINUE
  278. **    dnorma=dnorma*xzprec*xzprec
  279.       dnorma=xpetit/xzprec
  280.       SEGDES,MDIAG
  281.  
  282. C
  283. C  ****  MONTEE  ****
  284. C
  285. C ... Si la matrice est asymétrique ...
  286.       IF(EXISTS) THEN
  287.  
  288. C    ... On commence par désactiver les LIGN qui sont restés actifs ...
  289.          DO 2000 I=IDEP,ILMAX
  290.             LIGN=ILIGN(I)
  291.             if (i.lt.ilmin) seGDES,LIGN*(NOMOD,MRU)
  292.  2000    CONTINUE
  293.  
  294. C    ... Puis on désactive MILIGN ...
  295.          SEGDES,MILIGN
  296.  
  297. C    ... On change de MILIGN ...
  298.          MILIGN=IILIGS
  299.          SEGACT,MILIGN
  300.  
  301. C    ... Et on active des LIGN pointés par ce dernier ...
  302.          ILMAX=IDEP-1
  303.          DO 2100 I=IDEP,NNOE
  304.             LIGN=ILIGN(I)
  305.             SEGACT /ERR=2110/ LIGN
  306.             ILMAX = I
  307.  2100    CONTINUE
  308. C    ... Et on passe à la montée ...
  309.          GOTO 3000
  310.  
  311. C    ... Si on n'a même pas activé le premier, alors ADIOS ...
  312.  2110    IF(ILMAX.EQ.IDEP-1) CALL ERREUR(5)
  313.          toutmem=.false.
  314.  3000    CONTINUE
  315.       ENDIF
  316.  
  317.  
  318.       J=NNOE
  319. C ... Début de la pseudo boucle sur les LIGN qui sont restés désactivés ...
  320.   70  CONTINUE
  321.          IF (J.EQ.ILMAX) GOTO 72
  322.          LIGN=ILIGN(J)
  323.          SEGACT /ERR=68/ LIGN
  324.          GO TO 67
  325. C    ... Si on a des pb avec activation, on désactive et on diminue ILMAX ...
  326.  68      CONTINUE
  327.          toutmem=.false.
  328.          LIGN = ILIGN( ILMAX)
  329.          if(ilmax.lt.ilmin) SEGDES LIGN*(NOMOD,MRU)
  330.          ILMAX=ILMAX-1
  331.          IF(ILMAX.EQ.1) CALL ERREUR (5)
  332. C    ... puis on recommence la tentative ...
  333.          GO TO 70
  334. C    ... On a réussi à activer ...
  335.   67     CONTINUE
  336.          NA=IMMM(/1)
  337.          IPRELL=IPREL
  338.          CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  339.      >     NA,IPREL,MULRE,INC,dnorma)
  340.          if (j.lt.ilmin) SEGDES,LIGN*(NOMOD,MRU)
  341.          J = J-1
  342.          GO TO 70
  343.   72  CONTINUE
  344. CC  FIN DE PSEUDO BOUCLE J = INC ,ILMAX+1,-1 - Vieux commentaire et Faux !!!
  345. CC  FIN DE PSEUDO BOUCLE J = NNOE ,ILMAX+1,-1
  346.  
  347. C ... Dans cette boucle on parcourt les LIGN qui sont restés actifs ...
  348.       DO 13 J=ILMAX,IDEP,-1
  349.          LIGN=ILIGN(J)
  350.          NA=IMMM(/1)
  351.          IPRELL=IPREL
  352.          CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  353.      >     NA,IPREL,MULRE,INC,dnorma)
  354.          if (j.lt.ilmin) SEGDES,LIGN*(NOMOD,MRU)
  355.    13 CONTINUE
  356.  
  357. C ... On n'avait même pas essayé d'activer les IDEP-1 premiers LIGN ...
  358.       DO 113 J=IDEP-1,1,-1
  359.          LIGN=ILIGN(J)
  360.          SEGACT LIGN
  361.          NA=IMMM(/1)
  362.          IPRELL=IPREL
  363.          DO 1131 ILM=1,NA
  364.             IF(IMMM(ILM).EQ.0)  GO TO 1131
  365.             IAA=IAA+1
  366.             ITRAA(IAA)=IPRELL+ILM-1
  367.  1131    CONTINUE
  368.          CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  369.      >     NA,IPREL,MULRE,INC,dnorma)
  370.          if (j.lt.ilmin) SEGDES,LIGN*(NOMOD,MRU)
  371.   113 CONTINUE
  372.  
  373.          if (idenorm.eq.1) then
  374.           call erreur(1049)
  375.           do k=0,muinc,inc
  376.           do i=1,inc
  377.              vectbb(i+k)=0.d0
  378.           enddo
  379.           enddo
  380.          endif
  381. *  desactivation generale
  382.        do ivs=1,nnoe
  383.           lign=ilign(ivs)
  384.           segdes lign
  385.        enddo
  386.  
  387. *
  388. *    on verifie la precision du resultat si pas de noen ou pas de modes d'ensemble
  389. **    write(6,*) ' noen nens ',noen,nens
  390.       if (jlicre.ne.0) then
  391.       if (noen.ne.0.or.nens.eq.0) then
  392.       if(mulre.eq.1.and..not.exists) then
  393. *  U+dU
  394.        if (icorrec.eq.0) then
  395.        do i=1,inc
  396.          mvect2.vectbb(i)=vectbb(i)
  397.        enddo
  398.        else
  399.        do i=1,inc
  400.          mvect2.vectbb(i)=mvect2.vectbb(i)+vectbb(i)
  401.        enddo
  402.        endif
  403.        ilicre=jlicre
  404.        segact ilicre
  405.        ligcre=ligcrp
  406.        segact ligcre
  407. *  K*U
  408. *  sauver l'ancien residu
  409.        if (mvect4.ne.0)    segsup mvect4   
  410.        if (icorrec.ne.0)    mvect4 = mvect3
  411.        segini mvect3
  412.        call graco7(ilicre,mvect2,mvect3,inc,nbthr,0)
  413. *  F-K*U
  414.        crit=xpetit/xzprec
  415.        crit2 = 0.D0
  416.        do i=1,inc
  417.         crit=max(crit,abs(mvect1.vectbb(i)),abs(mvect3.vectbb(i)))
  418.         mvect3.vectbb(i)=mvect1.vectbb(i)-mvect3.vectbb(i)
  419. **      write(6,*) ' mvect1 mvect3 ',mvect1.vectbb(i),mvect3.vectbb(i)
  420.         crit2=max(crit2,abs(mvect3.vectbb(i)))
  421.        enddo
  422. *  calcul lambda delta u
  423.        xlbd = 1.
  424.        if (icorrec.ne.0) then
  425.        xd = 0.
  426.        xn = 0.
  427.        do i=1,inc
  428.          r34         =mvect3.vectbb(i)-mvect4.vectbb(i)
  429.          xn = xn + mvect4.vectbb(i) * r34
  430.          xd = xd + r34 * r34
  431.        enddo
  432.        xd = max(xd,abs(xn)*xpetit,xpetit)
  433.        xlbd = -xn / xd
  434. **     write(6,*) ' lambda ',xlbd 
  435.        if (xlbd.lt.-1..or.xlbd.gt.2.) xlbd = 0.5
  436.        endif
  437. * test
  438.        if (icorrec.le.8) crite = crit * sqrt(xzprec*xszpre)
  439.  
  440. * critere lache apres les premieres ite©rations. demande une reflexion approfondie.
  441.        if (icorrec.ge.13) crite = crit * 1d-4
  442.        if(icorrec.gt.12.and.nensl.ne.0) crite=xgrand*xzprec
  443. **      write(6,*) ' crite crit2',crite,crit2,icorrec
  444.        if (crit2.gt.crite) then
  445. * il faut iterer
  446.          do i2=1,inc
  447.            vectbb(i2)=mvect3.vectbb(i2)* xlbd  
  448.          enddo
  449. **       segsup mvect3
  450.          icorrec=icorrec+1
  451.          if (icorrec.lt.15.and.toutmem) goto 10000
  452.        endif
  453. *  verif correcte
  454.        segsup mvect1,mvect3
  455.        if (mvect4.ne.0)    segsup mvect4   
  456.        do i=1,inc
  457.           vectbb(i)=mvect2.vectbb(i)
  458.        enddo
  459.        segsup mvect3
  460.        reaerr(1)=crit2/crit
  461.        interr(1)=icorrec
  462.        interr(2)=nensl
  463. ***    isouci = 0
  464.        if (icorrec.ge.15.or.crit2/crit.gt.1d-7) then
  465.          if (isouci.eq.0) then
  466.             call erreur(1128)
  467.          else
  468.             call soucis(1128)
  469.             call erreur(1129)
  470.          endif
  471.        endif
  472.        if (icorrec.ge.1) then
  473. *  l'erreur 1129 est informative donc pas de souci
  474. **          call erreur(1129)
  475.        endif
  476.       endif
  477.       endif
  478.       endif
  479. **    write(6,*) 'icorrec',icorrec
  480. C
  481. C  ****  ON VERIFIE QUE PAS DE MODE RIGIDE ACTIVE
  482. C
  483. C ... IAA = nombre de modes d'ensemble (selon les indications dans IMMM) ...
  484.       NENS=IAA
  485. **    write(6,*) ' mondes nens nensl ',nens,nensl
  486.       NBEN=0
  487.  
  488.       MINCPO=IINCPO
  489.       MIMIK=IIMIK
  490.       SEGACT,MINCPO,MIMIK
  491.       MELEME=IGEOMA
  492.       SEGACT,MELEME
  493.       NNOE=INCPO(/2)
  494.       IINC1=INCPO(/1)
  495.  
  496. *  Pour creation chpoint des modes d'ensemble actifs
  497.       IF (NOEN.EQ.0) then
  498.        NNIN=IMIK(/2)
  499.        NNNOE=NUM(/2)
  500.        SEGINI MTRAV
  501.       DO I=1,IMIK(/2)
  502.        INCO(I)=IMIK(I)
  503.       ENDDO
  504.       ENDIF
  505. **     write(6,*) 'mondes creation mtrav:',mtrav,noen
  506.       IF(NENS.EQ.0)  GO TO 26
  507.  
  508. C ... Boucle sur les seconds membres ...
  509.       DO 300 KI=0,MUINC,INC
  510. C    ... XMA et XMAL seront le maximum des val. abs. des termes
  511. C        du résultat partiel (avant la multiplication par MDNOR)
  512. C        N° KI+1 multiplié par 1.e-10 ...
  513.          XMA=xpetit/xzprec
  514.          XMAL=xpetit/xzprec
  515.          DO 30 KK=1,INC
  516.             I=KK+KI
  517.             if (ittr(kk).eq.0) then
  518.              XMA=MAX(XMA,ABS(VECTBB(I)))
  519.             else
  520.              XMAL=MAX(XMAL,ABS(VECTBB(I)))
  521.             endif
  522.    30    CONTINUE
  523.          XMA=XMA*1.d-10
  524.          XMAL=XMAL*max(1d-10,crit2/crit*1d5)
  525.          xmal=max(xma*1d-2,xmal)
  526. *        write (6,*) ' mondes xma xmal ',xma,xmal
  527.  
  528. C    ... Boucle sur les modes d'ensemble ...
  529. **       write (6,*) ' nombre de modes d ensemble',nens
  530.          iwrite = 0
  531.          DO 21 IA=1,NENS
  532.             I1=ITRAA(IA)
  533.             J=IPNO(I1)
  534.             DO 22 K=1,IINC1
  535.                IF(INCPO(K,J).EQ.I1) GO TO 23
  536.    22       CONTINUE
  537.                CALL ERREUR(5)
  538.                call oooprl(0)
  539.                RETURN
  540.  
  541.    23       CONTINUE
  542. C       ... Si ce n'est pas un multiplicateur, le déplacement doit être
  543. C           < à XMA, sinon ERREUR ...
  544. *           write (6,*) ' mondes i1 ittr vect xma ',
  545. *    >             i1,ittr(i1),vectbb(i1+ki),xma
  546.             IF(ITTR(I1).EQ.0.AND.ABS(VECTBB(I1+KI)).LE.XMA    ) then
  547.               vectbb(i1+ki)=0.d0
  548.               GO TO 20
  549.             endif
  550. C       ... Si c'est un multiplicateur, le multiplicateur doit être
  551. C           < à XMAL sinon ERREUR ...
  552.             i2=1
  553.             if (ittr(i1).ne.0) then
  554.              i2=ittr(i1)
  555.              vecs=vectbb(i1+ki)+vectbb(i2+ki)
  556.              if (abs(vecs).le.xmal) then
  557.               vectbb(i1+ki)=0.d0
  558.               vectbb(i2+ki)=0.d0
  559.               goto 20
  560.              endif
  561.             endif
  562. *           write (6,*) ' ittr vect ',i1,i2,vectbb(i1+ki),
  563. *    >      vectbb(i2+ki),k,imik(k)
  564.  
  565.  
  566. C ... Si option NOEN on ne fait pas d'erreur, on accumule les pts dans un meleme
  567. C ... Si option NOEN on accumule les pts dans un meleme
  568.             IF(NOEN.EQ.0) THEN
  569.              IBIN(k,J)=1
  570.              IGEO(J)=num(1,j)
  571.              vv = vectbb(i1+ki)
  572.              bb(k,j)=vectbb(i1+ki)
  573. *  mettre aussi le second multiplicateur de lagrange pour l'appel a dbbcf
  574.              if (ittr(i1).ne.0) then
  575.              jd = ipno(i2)
  576.              IBIN(k,Jd)=1
  577.              IGEO(Jd)=num(1,jd)
  578.              vv = vectbb(i1+ki)
  579.              bb(k,jd)=vectbb(i2+ki)
  580.              endif
  581.             ENDIF
  582.                NBEN=NBEN+1
  583.                MOTERR(1:4)=IMIK(K)
  584.                INTERR(1)=NUM(1,J)
  585.                IF(NOEN.EQ.0.OR.ISOUCI.EQ.1) THEN
  586.                  call soucis(149)
  587.                  GO TO 21
  588.                ENDIF
  589.                CALL ERREUR(149)
  590.                call oooprl(0)
  591.                RETURN
  592.  
  593.    20       CONTINUE
  594. C       ... Messages d'information ...
  595.             JJK=NUM(1,J)
  596.             KIK=KI/INC +1
  597.  
  598.  
  599. *  on n'ecrit qu'une seule fois le message indetermination
  600.             IF(ITTR(I1).EQ.0 .AND. MULRE.EQ.1.and.iwrite.eq.0) then
  601.               if (imik(k).ne.'LX  ')
  602.      >                WRITE(IOIMP,41) JJK,IMIK(K)
  603.                 iwrite=iwrite+1
  604.             endif
  605.  
  606.             IF (IIMPI.NE.0.AND. ITTR(I1).NE.0.AND.MULRE.EQ.1) then
  607.                 WRITE(IOIMP,42) JJK,IMIK(K)
  608.             endif
  609.  
  610.             IF(ITTR(I1).eq.0 .AND. MULRE.NE.1.and.iwrite.eq.0) then
  611.                    WRITE(IOIMP,410)KIK,JJK,IMIK(K)
  612.                 iwrite=iwrite+1
  613.             endif
  614.             IF (IIMPI.NE.0 .AND. ITTR(I1).NE.0.AND.MULRE.NE.1)
  615.      &                WRITE(IOIMP,420)KIK,JJK,IMIK(K)
  616.  
  617.    21    CONTINUE
  618.  
  619.   300 CONTINUE
  620.  
  621.    41 FORMAT(' INDETERMINATION DETECTEE AU NOEUD ',I6,' INCONNUE ',
  622.      *   A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ',
  623.      * 'LA SUSDITE dans mondes')
  624.    42 FORMAT(' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ',
  625.      * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ',
  626.      * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE dans mondes')
  627.   410 FORMAT(' VECTEUR NUMERO',I3,' INDETERMINATION DETECTE AU NOEUD '
  628.      *,I6,' INCONNUE ',
  629.      *   A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ',
  630.      * 'LA SUSDITE dans mondes')
  631.   420 FORMAT(' VECTEUR NUMERO ',I3,/,
  632.      * ' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ',
  633.      * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ',
  634.      * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE dans mondes')
  635.  
  636.       SEGDES,MELEME
  637.       SEGDES,MINCPO
  638.       SEGDES,MIMIK
  639.    26 CONTINUE
  640.       if (noen.eq.0) then
  641. *  test si nan ou inf dans le resultat
  642.        inan = 0
  643.        DO 500 KI=0,MUINC,INC
  644.          DO 501 KK=1,INC
  645.           i = kk + ki
  646.           if (abs(vectbb(i)).lt.xgrand) goto 501
  647.           inan = inan + 1
  648.           vectbb(i)=xgrand
  649.  501   continue
  650.  500   continue
  651.        if (inan.ne.0) then
  652.          nben = -inan
  653.          call soucis(nben)
  654.        endif
  655.       endif
  656. *  indiquer si besoin le nombre de modes d'ensemble excites
  657. *  et le maillage des noeuds freinés
  658. *  et le champoint des modes excites
  659.       IF (NOEN.EQ.0) then
  660. *      if (nben.ne.0) write (6,*) ' mondes nben ',nben
  661.        CALL ECRENT(NBEN)
  662. *  pour permettre a crechp de prendre le persistent lock
  663.        call oooprl(0)
  664. **     write(6,*) 'mondes appeler crechp avec mtrav',mtrav
  665.        CALL CRECHP(MTRAV,MCHPO4)
  666.        call oooprl(1)
  667.        if(mtrav.ne.0) SEGSUP MTRAV
  668. *  dedualisation des multiplicateurs
  669.        if (lagdua.ne.0) call dbbcf(mchpo4,lagdua)
  670.        call ecrobj('CHPOINT ',mchpo4)
  671.       endif
  672.       MDNOR=IDNORM
  673.       SEGACT,MDNOR
  674.       DO 350 K=0,MUINC,INC
  675.          DO 35 I = 1, INC
  676.             VECTBB(I+K)=VECTBB(I+K)*DNOR(I)
  677.    35    CONTINUE
  678. *  on force l'egalite des multiplicateurs de lagrange
  679.          DO 36 I = 1, INC
  680.             if (ITTR(I).ne.0) then
  681.              vectbb(i+k)=(vectbb(i+k)+vectbb(ittr(i)+k))*0.5D0
  682.              vectbb(ittr(i)+k)=vectbb(i+k)
  683.             endif
  684.    36    CONTINUE
  685.   350 CONTINUE
  686.       SEGDES,MDNOR
  687. C ... On ne désactive MDNO1 que dans le cas où il est
  688. C     différent de MDNOR ...
  689.       IF(IDNORD.GT.0) THEN
  690.          SEGDES,MDNO1
  691.       ENDIF
  692.       SEGDES MMATRI
  693.       SEGDES,MILIGN
  694.       SEGDES,MVECTD
  695.       SEGSUP ITTRV
  696.       IF(NENSL.NE.0) SEGSUP,ITRAA
  697.  
  698.       IF(IIMPI.EQ.3) WRITE(IOIMP,1001) MVECTD
  699.  1001 FORMAT('  SUBROUTINE MONDES : POINTEUR DU VECTEUR EN SORTIE=',I5)
  700.  
  701.       CALL OOOMRU(0)
  702.       call oooprl(0)
  703.       RETURN
  704.       END
  705.  
  706.  
  707.  
  708.  
  709.  
  710.  
  711.  
  712.  
  713.  
  714.  
  715.  
  716.  
  717.  
  718.  
  719.  
  720.  
  721.  
  722.  
  723.  
  724.  
  725.  
  726.  
  727.  
  728.  
  729.  
  730.  
  731.  
  732.  
  733.  
  734.  
  735.  
  736.  
  737.  
  738.  

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