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

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