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/07/24    14:05:24     12332          
  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=ilmin,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 (noen.ne.0.or.nens.eq.0) then
  391.       if(mulre.eq.1.and..not.exists) then
  392. *  U+dU
  393.        if (icorrec.eq.0) then
  394.        do i=1,inc
  395.          mvect2.vectbb(i)=vectbb(i)
  396.        enddo
  397.        else
  398.        do i=1,inc
  399.          mvect2.vectbb(i)=mvect2.vectbb(i)+vectbb(i)
  400.        enddo
  401.        endif
  402.        ilicre=jlicre
  403.        segact ilicre
  404.        ligcre=ligcrp
  405.        segact ligcre
  406. *  K*U
  407. *  sauver l'ancien residu
  408.        if (mvect4.ne.0)    segsup mvect4   
  409.        if (icorrec.ne.0)    mvect4 = mvect3
  410.        segini mvect3
  411.        call graco7(ilicre,mvect2,mvect3,inc,nbthr,0)
  412. *  F-K*U
  413.        crit=xpetit/xzprec
  414.        crit2 = 0.D0
  415.        do i=1,inc
  416.         crit=max(crit,abs(mvect1.vectbb(i)),abs(mvect3.vectbb(i)))
  417.         mvect3.vectbb(i)=mvect1.vectbb(i)-mvect3.vectbb(i)
  418. **      write(6,*) ' mvect1 mvect3 ',mvect1.vectbb(i),mvect3.vectbb(i)
  419.         crit2=max(crit2,abs(mvect3.vectbb(i)))
  420.        enddo
  421. *  calcul lambda delta u
  422.        xlbd = 1.
  423.        if (icorrec.ne.0) then
  424.        xd = 0.
  425.        xn = 0.
  426.        do i=1,inc
  427.          r34         =mvect3.vectbb(i)-mvect4.vectbb(i)
  428.          xn = xn + mvect4.vectbb(i) * r34
  429.          xd = xd + r34 * r34
  430.        enddo
  431.        xd = max(xd,abs(xn)*xpetit,xpetit)
  432.        xlbd = -xn / xd
  433. **     write(6,*) ' lambda ',xlbd 
  434.        if (xlbd.lt.-1..or.xlbd.gt.2.) xlbd = 0.5
  435.        endif
  436. * test
  437.        if (icorrec.le.8) crite = crit * sqrt(xzprec*xszpre)
  438.  
  439. * critere lache apres les premieres ite©rations. demande une reflexion approfondie.
  440.        if (icorrec.ge.13) crite = crit * 1d-4
  441.        if(icorrec.gt.12.and.nensl.ne.0) crite=xgrand*xzprec
  442. **      write(6,*) ' crite crit2',crite,crit2,icorrec
  443.        if (crit2.gt.crite) then
  444. * il faut iterer
  445.          do i2=1,inc
  446.            vectbb(i2)=mvect3.vectbb(i2)* xlbd  
  447.          enddo
  448. **       segsup mvect3
  449.          icorrec=icorrec+1
  450.          if (icorrec.lt.15.and.toutmem) goto 10000
  451.        endif
  452. *  verif correcte
  453.        segsup mvect1,mvect3
  454.        if (mvect4.ne.0)    segsup mvect4   
  455.        do i=1,inc
  456.           vectbb(i)=mvect2.vectbb(i)
  457.        enddo
  458.        segsup mvect3
  459.        reaerr(1)=crit2/crit
  460.        interr(1)=icorrec
  461.        interr(2)=nensl
  462. ***    isouci = 0
  463.        if (icorrec.ge.15.or.crit2/crit.gt.1d-7) then
  464.          if (isouci.eq.0) then
  465.             call erreur(1128)
  466.          else
  467.             call soucis(1128)
  468.             call erreur(1129)
  469.          endif
  470.        endif
  471.        if (icorrec.ge.1) then
  472. *  l'erreur 1129 est informative donc pas de souci
  473. **          call erreur(1129)
  474.        endif
  475.       endif
  476.       endif
  477. **    write(6,*) 'icorrec',icorrec
  478. C
  479. C  ****  ON VERIFIE QUE PAS DE MODE RIGIDE ACTIVE
  480. C
  481. C ... IAA = nombre de modes d'ensemble (selon les indications dans IMMM) ...
  482.       NENS=IAA
  483. **    write(6,*) ' mondes nens nensl ',nens,nensl
  484.       NBEN=0
  485.  
  486.       MINCPO=IINCPO
  487.       MIMIK=IIMIK
  488.       SEGACT,MINCPO,MIMIK
  489.       MELEME=IGEOMA
  490.       SEGACT,MELEME
  491.       NNOE=INCPO(/2)
  492.       IINC1=INCPO(/1)
  493.  
  494. *  Pour creation chpoint des modes d'ensemble actifs
  495.       IF (NOEN.EQ.0) then
  496.        NNIN=IMIK(/2)
  497.        NNNOE=NUM(/2)
  498.        SEGINI MTRAV
  499.       DO I=1,IMIK(/2)
  500.        INCO(I)=IMIK(I)
  501.       ENDDO
  502.       ENDIF
  503. **     write(6,*) 'mondes creation mtrav:',mtrav,noen
  504.       IF(NENS.EQ.0)  GO TO 26
  505.  
  506. C ... Boucle sur les seconds membres ...
  507.       DO 300 KI=0,MUINC,INC
  508. C    ... XMA et XMAL seront le maximum des val. abs. des termes
  509. C        du résultat partiel (avant la multiplication par MDNOR)
  510. C        N° KI+1 multiplié par 1.e-10 ...
  511.          XMA=xpetit/xzprec
  512.          XMAL=xpetit/xzprec
  513.          DO 30 KK=1,INC
  514.             I=KK+KI
  515.             if (ittr(kk).eq.0) then
  516.              XMA=MAX(XMA,ABS(VECTBB(I)))
  517.             else
  518.              XMAL=MAX(XMAL,ABS(VECTBB(I)))
  519.             endif
  520.    30    CONTINUE
  521.          XMA=XMA*1.d-10
  522.          XMAL=XMAL*max(1d-10,crit2/crit*1d5)
  523.          xmal=max(xma*1d-2,xmal)
  524. *        write (6,*) ' mondes xma xmal ',xma,xmal
  525.  
  526. C    ... Boucle sur les modes d'ensemble ...
  527. **       write (6,*) ' nombre de modes d ensemble',nens
  528.          iwrite = 0
  529.          DO 21 IA=1,NENS
  530.             I1=ITRAA(IA)
  531.             J=IPNO(I1)
  532.             DO 22 K=1,IINC1
  533.                IF(INCPO(K,J).EQ.I1) GO TO 23
  534.    22       CONTINUE
  535.                CALL ERREUR(5)
  536.                call oooprl(0)
  537.                RETURN
  538.  
  539.    23       CONTINUE
  540. C       ... Si ce n'est pas un multiplicateur, le déplacement doit être
  541. C           < à XMA, sinon ERREUR ...
  542. *           write (6,*) ' mondes i1 ittr vect xma ',
  543. *    >             i1,ittr(i1),vectbb(i1+ki),xma
  544.             IF(ITTR(I1).EQ.0.AND.ABS(VECTBB(I1+KI)).LE.XMA    ) then
  545.               vectbb(i1+ki)=0.d0
  546.               GO TO 20
  547.             endif
  548. C       ... Si c'est un multiplicateur, le multiplicateur doit être
  549. C           < à XMAL sinon ERREUR ...
  550.             i2=1
  551.             if (ittr(i1).ne.0) then
  552.              i2=ittr(i1)
  553.              vecs=vectbb(i1+ki)+vectbb(i2+ki)
  554.              if (abs(vecs).le.xmal) then
  555.               vectbb(i1+ki)=0.d0
  556.               vectbb(i2+ki)=0.d0
  557.               goto 20
  558.              endif
  559.             endif
  560. *           write (6,*) ' ittr vect ',i1,i2,vectbb(i1+ki),
  561. *    >      vectbb(i2+ki),k,imik(k)
  562.  
  563.  
  564. C ... Si option NOEN on ne fait pas d'erreur, on accumule les pts dans un meleme
  565. C ... Si option NOEN on accumule les pts dans un meleme
  566.             IF(NOEN.EQ.0) THEN
  567.              IBIN(k,J)=1
  568.              IGEO(J)=num(1,j)
  569.              vv = vectbb(i1+ki)
  570.              bb(k,j)=vectbb(i1+ki)
  571. *  mettre aussi le second multiplicateur de lagrange pour l'appel a dbbcf
  572.              if (ittr(i1).ne.0) then
  573.              jd = ipno(i2)
  574.              IBIN(k,Jd)=1
  575.              IGEO(Jd)=num(1,jd)
  576.              vv = vectbb(i1+ki)
  577.              bb(k,jd)=vectbb(i2+ki)
  578.              endif
  579.             ENDIF
  580.                NBEN=NBEN+1
  581.                MOTERR(1:4)=IMIK(K)
  582.                INTERR(1)=NUM(1,J)
  583.                IF(NOEN.EQ.0.OR.ISOUCI.EQ.1) THEN
  584.                  call soucis(149)
  585.                  GO TO 21
  586.                ENDIF
  587.                CALL ERREUR(149)
  588.                call oooprl(0)
  589.                RETURN
  590.  
  591.    20       CONTINUE
  592. C       ... Messages d'information ...
  593.             JJK=NUM(1,J)
  594.             KIK=KI/INC +1
  595.  
  596.  
  597. *  on n'ecrit qu'une seule fois le message indetermination
  598.             IF(ITTR(I1).EQ.0 .AND. MULRE.EQ.1.and.iwrite.eq.0) then
  599.               if (imik(k).ne.'LX  ')
  600.      >                WRITE(IOIMP,41) JJK,IMIK(K)
  601.                 iwrite=iwrite+1
  602.             endif
  603.  
  604.             IF (IIMPI.NE.0.AND. ITTR(I1).NE.0.AND.MULRE.EQ.1) then
  605.                 WRITE(IOIMP,42) JJK,IMIK(K)
  606.             endif
  607.  
  608.             IF(ITTR(I1).eq.0 .AND. MULRE.NE.1.and.iwrite.eq.0) then
  609.                    WRITE(IOIMP,410)KIK,JJK,IMIK(K)
  610.                 iwrite=iwrite+1
  611.             endif
  612.             IF (IIMPI.NE.0 .AND. ITTR(I1).NE.0.AND.MULRE.NE.1)
  613.      &                WRITE(IOIMP,420)KIK,JJK,IMIK(K)
  614.  
  615.    21    CONTINUE
  616.  
  617.   300 CONTINUE
  618.  
  619.    41 FORMAT(' INDETERMINATION DETECTEE AU NOEUD ',I6,' INCONNUE ',
  620.      *   A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ',
  621.      * 'LA SUSDITE dans mondes')
  622.    42 FORMAT(' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ',
  623.      * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ',
  624.      * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE dans mondes')
  625.   410 FORMAT(' VECTEUR NUMERO',I3,' INDETERMINATION DETECTE AU NOEUD '
  626.      *,I6,' INCONNUE ',
  627.      *   A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ',
  628.      * 'LA SUSDITE dans mondes')
  629.   420 FORMAT(' VECTEUR NUMERO ',I3,/,
  630.      * ' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ',
  631.      * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ',
  632.      * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE dans mondes')
  633.  
  634.       SEGDES,MELEME
  635.       SEGDES,MINCPO
  636.       SEGDES,MIMIK
  637.    26 CONTINUE
  638.       if (noen.eq.0) then
  639. *  test si nan ou inf dans le resultat
  640.        inan = 0
  641.        DO 500 KI=0,MUINC,INC
  642.          DO 501 KK=1,INC
  643.           i = kk + ki
  644.           if (abs(vectbb(i)).lt.xgrand) goto 501
  645.           inan = inan + 1
  646.           vectbb(i)=xgrand
  647.  501   continue
  648.  500   continue
  649.        if (inan.ne.0) then
  650.          nben = -inan
  651.          call soucis(nben)
  652.        endif
  653.       endif
  654. *  indiquer si besoin le nombre de modes d'ensemble excites
  655. *  et le maillage des noeuds freinés
  656. *  et le champoint des modes excites
  657.       IF (NOEN.EQ.0) then
  658. *      if (nben.ne.0) write (6,*) ' mondes nben ',nben
  659.        CALL ECRENT(NBEN)
  660. *  pour permettre a crechp de prendre le persistent lock
  661.        call oooprl(0)
  662. **     write(6,*) 'mondes appeler crechp avec mtrav',mtrav
  663.        CALL CRECHP(MTRAV,MCHPO4)
  664.        call oooprl(1)
  665.        if(mtrav.ne.0) SEGSUP MTRAV
  666. *  dedualisation des multiplicateurs
  667.        if (lagdua.ne.0) call dbbcf(mchpo4,lagdua)
  668.        call ecrobj('CHPOINT ',mchpo4)
  669.       endif
  670.       MDNOR=IDNORM
  671.       SEGACT,MDNOR
  672.       DO 350 K=0,MUINC,INC
  673.          DO 35 I = 1, INC
  674.             VECTBB(I+K)=VECTBB(I+K)*DNOR(I)
  675.    35    CONTINUE
  676. *  on force l'egalite des multiplicateurs de lagrange
  677.          DO 36 I = 1, INC
  678.             if (ITTR(I).ne.0) then
  679.              vectbb(i+k)=(vectbb(i+k)+vectbb(ittr(i)+k))*0.5D0
  680.              vectbb(ittr(i)+k)=vectbb(i+k)
  681.             endif
  682.    36    CONTINUE
  683.   350 CONTINUE
  684.       SEGDES,MDNOR
  685. C ... On ne désactive MDNO1 que dans le cas où il est
  686. C     différent de MDNOR ...
  687.       IF(IDNORD.GT.0) THEN
  688.          SEGDES,MDNO1
  689.       ENDIF
  690.       SEGDES MMATRI
  691.       SEGDES,MILIGN
  692.       SEGDES,MVECTD
  693.       SEGSUP ITTRV
  694.       IF(NENSL.NE.0) SEGSUP,ITRAA
  695.  
  696.       IF(IIMPI.EQ.3) WRITE(IOIMP,1001) MVECTD
  697.  1001 FORMAT('  SUBROUTINE MONDES : POINTEUR DU VECTEUR EN SORTIE=',I5)
  698.  
  699.       CALL OOOMRU(0)
  700.       call oooprl(0)
  701.       RETURN
  702.       END
  703.  
  704.  
  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.  

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