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  25/01/15    21:15:04     12125          
  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.  
  464.       MINCPO=IINCPO
  465.       MIMIK=IIMIK
  466.       SEGACT,MINCPO,MIMIK
  467.       MELEME=IGEOMA
  468.       SEGACT,MELEME
  469.       NNOE=INCPO(/2)
  470.       IINC1=INCPO(/1)
  471.  
  472. *  Pour creation chpoint des modes d'ensemble actifs
  473.       IF (NOEN.EQ.0) then
  474.        NNIN=IMIK(/2)
  475.        NNNOE=NUM(/2)
  476.        SEGINI MTRAV
  477.       DO I=1,IMIK(/2)
  478.        INCO(I)=IMIK(I)
  479.       ENDDO
  480.       ENDIF
  481. **     write(6,*) 'mondes creation mtrav:',mtrav,noen
  482.       IF(NENS.EQ.0)  GO TO 26
  483.  
  484. C ... Boucle sur les seconds membres ...
  485.       DO 300 KI=0,MUINC,INC
  486. C    ... XMA et XMAL seront le maximum des val. abs. des termes
  487. C        du résultat partiel (avant la multiplication par MDNOR)
  488. C        N° KI+1 multiplié par 1.e-10 ...
  489.          XMA=xpetit/xzprec
  490.          XMAL=xpetit/xzprec
  491.          DO 30 KK=1,INC
  492.             I=KK+KI
  493.             if (ittr(kk).eq.0) then
  494.              XMA=MAX(XMA,ABS(VECTBB(I)))
  495.             else
  496.              XMAL=MAX(XMAL,ABS(VECTBB(I)))
  497.             endif
  498.    30    CONTINUE
  499.          XMA=XMA*1.d-10
  500.          XMAL=XMAL*max(1d-10,crit2/crit*1d5)
  501.          xmal=max(xma*1d-2,xmal)
  502. *        write (6,*) ' mondes xma xmal ',xma,xmal
  503.  
  504. C    ... Boucle sur les modes d'ensemble ...
  505. **       write (6,*) ' nombre de modes d ensemble',nens
  506.          iwrite = 0
  507.          DO 21 IA=1,NENS
  508.             I1=ITRAA(IA)
  509.             J=IPNO(I1)
  510.             DO 22 K=1,IINC1
  511.                IF(INCPO(K,J).EQ.I1) GO TO 23
  512.    22       CONTINUE
  513.                CALL ERREUR(5)
  514.                call oooprl(0)
  515.                RETURN
  516.  
  517.    23       CONTINUE
  518. C       ... Si ce n'est pas un multiplicateur, le déplacement doit être
  519. C           < à XMA, sinon ERREUR ...
  520. *           write (6,*) ' mondes i1 ittr vect xma ',
  521. *    >             i1,ittr(i1),vectbb(i1+ki),xma
  522.             IF(ITTR(I1).EQ.0.AND.ABS(VECTBB(I1+KI)).LE.XMA    ) then
  523.               vectbb(i1+ki)=0.d0
  524.               GO TO 20
  525.             endif
  526. C       ... Si c'est un multiplicateur, le multiplicateur doit être
  527. C           < à XMAL sinon ERREUR ...
  528.             i2=1
  529.             if (ittr(i1).ne.0) then
  530.              i2=ittr(i1)
  531.              vecs=vectbb(i1+ki)+vectbb(i2+ki)
  532.              if (abs(vecs).le.xmal) then
  533.               vectbb(i1+ki)=0.d0
  534.               vectbb(i2+ki)=0.d0
  535.               goto 20
  536.              endif
  537.             endif
  538. *           write (6,*) ' ittr vect ',i1,i2,vectbb(i1+ki),
  539. *    >      vectbb(i2+ki),k,imik(k)
  540.  
  541.  
  542. C ... Si option NOEN on ne fait pas d'erreur, on accumule les pts dans un meleme
  543. C ... Si option NOEN on accumule les pts dans un meleme
  544.             IF(NOEN.EQ.0) THEN
  545.              IBIN(k,J)=1
  546.              IGEO(J)=num(1,j)
  547.              vv = vectbb(i1+ki)
  548.              bb(k,j)=vectbb(i1+ki)
  549. *  mettre aussi le second multiplicateur de lagrange pour l'appel a dbbcf
  550.              if (ittr(i1).ne.0) then
  551.              jd = ipno(i2)
  552.              IBIN(k,Jd)=1
  553.              IGEO(Jd)=num(1,jd)
  554.              vv = vectbb(i1+ki)
  555.              bb(k,jd)=vectbb(i2+ki)
  556.              endif
  557.             ENDIF
  558.                NBEN=NBEN+1
  559.                MOTERR(1:4)=IMIK(K)
  560.                INTERR(1)=NUM(1,J)
  561.                IF(NOEN.EQ.0.OR.ISOUCI.EQ.1) THEN
  562.                  call soucis(149)
  563.                  GO TO 21
  564.                ENDIF
  565.                CALL ERREUR(149)
  566.                call oooprl(0)
  567.                RETURN
  568.  
  569.    20       CONTINUE
  570. C       ... Messages d'information ...
  571.             JJK=NUM(1,J)
  572.             KIK=KI/INC +1
  573.  
  574.  
  575. *  on n'ecrit qu'une seule fois le message indetermination
  576.             IF(ITTR(I1).EQ.0 .AND. MULRE.EQ.1.and.iwrite.eq.0) then
  577.               if (imik(k).ne.'LX  ')
  578.      >                WRITE(IOIMP,41) JJK,IMIK(K)
  579.                 iwrite=iwrite+1
  580.             endif
  581.  
  582.             IF (IIMPI.NE.0.AND. ITTR(I1).NE.0.AND.MULRE.EQ.1) then
  583.                 WRITE(IOIMP,42) JJK,IMIK(K)
  584.             endif
  585.  
  586.             IF(ITTR(I1).eq.0 .AND. MULRE.NE.1.and.iwrite.eq.0) then
  587.                    WRITE(IOIMP,410)KIK,JJK,IMIK(K)
  588.                 iwrite=iwrite+1
  589.             endif
  590.             IF (IIMPI.NE.0 .AND. ITTR(I1).NE.0.AND.MULRE.NE.1)
  591.      &                WRITE(IOIMP,420)KIK,JJK,IMIK(K)
  592.  
  593.    21    CONTINUE
  594.  
  595.   300 CONTINUE
  596.  
  597.    41 FORMAT(' INDETERMINATION DETECTEE AU NOEUD ',I6,' INCONNUE ',
  598.      *   A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ',
  599.      * 'LA SUSDITE dans mondes')
  600.    42 FORMAT(' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ',
  601.      * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ',
  602.      * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE dans mondes')
  603.   410 FORMAT(' VECTEUR NUMERO',I3,' INDETERMINATION DETECTE AU NOEUD '
  604.      *,I6,' INCONNUE ',
  605.      *   A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ',
  606.      * 'LA SUSDITE dans mondes')
  607.   420 FORMAT(' VECTEUR NUMERO ',I3,/,
  608.      * ' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ',
  609.      * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ',
  610.      * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE dans mondes')
  611.  
  612.       SEGDES,MELEME
  613.       SEGDES,MINCPO
  614.       SEGDES,MIMIK
  615.    26 CONTINUE
  616.       if (noen.eq.0) then
  617. *  test si nan ou inf dans le resultat
  618.        inan = 0
  619.        DO 500 KI=0,MUINC,INC
  620.          DO 501 KK=1,INC
  621.           i = kk + ki
  622.           if (abs(vectbb(i)).lt.xgrand) goto 501
  623.           inan = inan + 1
  624.           vectbb(i)=xgrand
  625.  501   continue
  626.  500   continue
  627.        if (inan.ne.0) then
  628.          nben = -inan
  629.          call soucis(nben)
  630.        endif
  631.       endif
  632. *  indiquer si besoin le nombre de modes d'ensemble excites
  633. *  et le maillage des noeuds freinés
  634. *  et le champoint des modes excites
  635.       IF (NOEN.EQ.0) then
  636. *      if (nben.ne.0) write (6,*) ' mondes nben ',nben
  637.        CALL ECRENT(NBEN)
  638. *  pour permettre a crechp de prendre le persistent lock
  639.        call oooprl(0)
  640. **     write(6,*) 'mondes appeler crechp avec mtrav',mtrav
  641.        CALL CRECHP(MTRAV,MCHPO4)
  642.        call oooprl(1)
  643.        if(mtrav.ne.0) SEGSUP MTRAV
  644. *  dedualisation des multiplicateurs
  645.        if (lagdua.ne.0) call dbbcf(mchpo4,lagdua)
  646.        call ecrobj('CHPOINT ',mchpo4)
  647.       endif
  648.       MDNOR=IDNORM
  649.       SEGACT,MDNOR
  650.       DO 350 K=0,MUINC,INC
  651.          DO 35 I = 1, INC
  652.             VECTBB(I+K)=VECTBB(I+K)*DNOR(I)
  653.    35    CONTINUE
  654. *  on force l'egalite des multiplicateurs de lagrange
  655.          DO 36 I = 1, INC
  656.             if (ITTR(I).ne.0) then
  657.              vectbb(i+k)=(vectbb(i+k)+vectbb(ittr(i)+k))*0.5D0
  658.              vectbb(ittr(i)+k)=vectbb(i+k)
  659.             endif
  660.    36    CONTINUE
  661.   350 CONTINUE
  662.       SEGDES,MDNOR
  663. C ... On ne désactive MDNO1 que dans le cas où il est
  664. C     différent de MDNOR ...
  665.       IF(IDNORD.GT.0) THEN
  666.          SEGDES,MDNO1
  667.       ENDIF
  668.       SEGDES MMATRI
  669.       SEGDES,MILIGN
  670.       SEGDES,MVECTD
  671.       SEGSUP ITTRV
  672.       IF(NENSL.NE.0) SEGSUP,ITRAA
  673.  
  674.       IF(IIMPI.EQ.3) WRITE(IOIMP,1001) MVECTD
  675.  1001 FORMAT('  SUBROUTINE MONDES : POINTEUR DU VECTEUR EN SORTIE=',I5)
  676.  
  677.       CALL OOOMRU(0)
  678.       call oooprl(0)
  679.       RETURN
  680.       END
  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.  
  708.  
  709.  
  710.  

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