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

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