Télécharger la plaquette Cast3M


fr en

Télécharger excit1.eso

Retour à la liste

Numérotation des lignes : 

excit1

  1. C EXCIT1    SOURCE    MB234859  23/12/08    21:15:03     11806          
  2.  
  3.       SUBROUTINE EXCIT1(MRIGID,MCHPO1,MCHPO2,MCHPO3,MLENTI,indic,RI2,
  4.      &                  ipt8,ITYP)
  5.  
  6.       IMPLICIT INTEGER(I-N)
  7.       IMPLICIT REAL*8(A-H,O-Z)
  8.  
  9. *  modif : declarer inactifs les blocages redondants s'ils sont dans
  10. *  ipt8. Si celui-ci est non vide, on enleve les blocages inclus
  11. *  dedans et on s'en tient la.
  12.  
  13.  
  14. -INC PPARAM
  15. -INC CCOPTIO
  16. -INC CCREEL
  17. -INC CCGEOME
  18.  
  19. -INC SMRIGID
  20. -INC SMCHPOI
  21. -INC SMELEME
  22. -INC SMCOORD
  23. -INC SMLENTI
  24.       save icntr,icnts
  25.       DATA ICNTR/0/
  26.       DATA ICNTS/0/
  27.       SEGMENT SNOMIN
  28.         CHARACTER*(LOCOMP) NOMIN(0)
  29.       ENDSEGMENT
  30.       SEGMENT ISOU(IRIGEL(/2))
  31.       SEGMENT ICPR(nbpts)
  32.       SEGMENT JCPR(nbpts)
  33.       SEGMENT KCPR(nbpts)
  34.       SEGMENT YCPR(nbpts)
  35.       SEGMENT ITRAV
  36.         real*8 DVAL(NIN,NPO),DIN(NPO)
  37.         integer IEXI(NIN,NPO)
  38.       endsegment
  39.       SEGMENT  IPASS(ITEMP)
  40.       segment trav1
  41.         integer iav(nbo)
  42.         real*8 viol(nbo)
  43.       endsegment
  44. C
  45. C  ***  ICPR REFERENCIE DVAL
  46. C  ***  NOMIN CONTIENT LES INCONNUES REFERENCES PAR LES RELATIONS
  47. C  ***  DVAL CONTIENT LE RESULTAT DES DEPLACEMENTS ET DES LAMBDA
  48. C  ***  DIN CONTIENT LES VALEURS INITIALES DU SECOND MEMBRE
  49. C
  50.       logical indic,iclred
  51. C
  52.       idimp1 = idim+1
  53.       segact,mcoord
  54. C
  55.       XMCRIT = xpetit/xzprec
  56.       YMCRIT = xpetit/xzprec
  57.  
  58.       icpr=0
  59.       jcpr=0
  60.       ycpr=0
  61.       kcpr=0
  62.       itrav=0
  63.  
  64. * on a un maillage de relation(s) a enlever
  65.       if (ipt8.ne.0) then
  66.         segact,ipt8
  67.         if (ipt8.itypel.ne.1) then
  68.           call erreur(16)
  69.           return
  70.         endif
  71.         if (ipt8.num(/2).ne.0) then
  72.           segini,kcpr
  73.           do i=1,ipt8.num(/2)
  74.             kcpr(ipt8.num(1,i))=1
  75.           enddo
  76.         endif
  77.       endif
  78. *
  79. *  mchpo3 est le champ des LX sorti de excfro jcpr donne existence
  80. *   et ycpr donne valeur
  81. *     write (ioimp,*) ' mchpo3 dans excit1 ',mchpo3
  82.       if (mchpo3.ne.0) then
  83.         segini,jcpr,ycpr
  84.         segact,mchpo3
  85.         do 100 isoupo=1,mchpo3.ipchp(/1)
  86.           msoupo=mchpo3.ipchp(isoupo)
  87.           segact,msoupo
  88.           do 110 ic=1,nocomp(/2)
  89.             if (nocomp(ic).eq.'LX  ') goto 120
  90.  110      continue
  91.           goto 140
  92.  120      continue
  93.           ipt2=igeoc
  94.           mpoval=ipoval
  95.           segact,ipt2,mpoval
  96.           do 130 iel=1,ipt2.num(/2)
  97.             impt=ipt2.num(1,iel)
  98.             jcpr(impt)=1
  99.             ycpr(impt)=ycpr(impt)+vpocha(iel,ic)
  100.  130      continue
  101.  140      continue
  102.  100    continue
  103.       endif
  104. C
  105. C** INITIALISATION DE NOMIN ET ICPR icpr repere les points appartenant
  106. C** a la matrice de blocage mini maxi ou ?
  107. C
  108.       SEGINI,SNOMIN,ICPR
  109.       NOMIN(**)='LX  '
  110.       NPO=0
  111.       NBO=0
  112. *
  113.       SEGACT,MRIGID
  114.       DO 1 I=1,IRIGEL(/2)
  115.         ITY=IRIGEL(6,I)
  116.         IF (ITY.EQ.0) GO TO 1
  117. *  cas du frottement : petite verification
  118.         if (ity.eq.2 .and. mchpo3.eq.0) then
  119.           call erreur(721)
  120.           return
  121.         endif
  122.         MELEME=IRIGEL(1,I)
  123.         SEGACT,MELEME
  124.         NBO=NBO+NUM(/2)
  125.         IF(IIMPI.EQ.528) WRITE(IOIMP,3765) NBO
  126.  3765   FORMAT('  NBO  ',I5)
  127.         DO  J=1,NUM(/2)
  128.           DO  K=1,NUM(/1)
  129.             impt=NUM(K,J)
  130.             IF (ICPR(impt).EQ.0) THEN
  131.               NPO=NPO+1
  132.               ICPR(impt)=NPO
  133.             ENDIF
  134.           enddo
  135.         enddo
  136.         DESCR=IRIGEL(3,I)
  137.         SEGACT,DESCR
  138.         DO 3 J=1,LISINC(/2)
  139.           DO 4 K=1,NOMIN(/2)
  140.             IF (NOMIN(K).EQ.LISINC(J)) GO TO 3
  141.     4     CONTINUE
  142.           NOMIN(**)=LISINC(J)
  143.     3   CONTINUE
  144.         SEGDES,DESCR
  145.  1    CONTINUE
  146.       NIN=NOMIN(/2)
  147. *
  148.       IF(IIMPI.EQ.528) THEN
  149.         WRITE(ioimp,90) NPO
  150.         WRITE(IOIMP,9876)( NOMIN(i),i=1,NIN)
  151.   90    FORMAT(' ON VIENT DE PASSER LA BOUCLE 1  NPO  ',I5)
  152.  9876   FORMAT('  NOMIN ' ,10(A8,1X))
  153.       ENDIF
  154. C
  155. C ****  ON REMPLIT DVAL AVEC LES VALEURS DE MCHPO2
  156. C   mchpo2 est le champ de deplacement propose par RESOU
  157. C   dval(j,icpr(i)) est la j eme composante du deplacement du point i
  158. C   et iexi(j,icpr(i))=1
  159. C
  160.       SEGINI,ITRAV
  161. c
  162.       MCHPOI=MCHPO2
  163.       SEGACT,MCHPOI
  164.       DO 5 I=1,IPCHP(/1)
  165.         MSOUPO=IPCHP(I)
  166.         SEGACT,MSOUPO
  167.         MELEME=IGEOC
  168.         SEGACT,MELEME
  169.         ITEMP=NOCOMP(/2)
  170.         SEGINI,IPASS
  171.         DO 6 K=1,ITEMP
  172.           DO 7 J=1,NIN
  173.             IF (NOMIN(J).EQ.NOCOMP(K)) THEN
  174.               IPASS(K)=J
  175.               GO TO 6
  176.             ENDIF
  177.    7      CONTINUE
  178.    6    CONTINUE
  179.         IF (IIMPI.EQ.528) WRITE(IOIMP,1555) (IPASS(KHU),KHU=1,ITEMP)
  180.  1555 FORMAT('   IPASS   ' ,9I10)
  181.         MPOVAL=IPOVAL
  182.         SEGACT,MPOVAL
  183.         III=0
  184.         DO 8 J=1,ITEMP
  185.           K=IPASS(J)
  186.           IF (K.EQ.0) GO TO 8
  187.           DO 9 L=1,NUM(/2)
  188.             IP=ICPR(NUM(1,L))
  189.             IF (IP.EQ.0) GO TO 9
  190. *     if (k.eq.1.and.vpocha(l,j).eq.0.d0 ) then
  191. *      write (6,*) 'LX nul dans excit1'
  192. *      goto 9
  193. *     endif
  194.             DVAL(K,IP)=dval(k,ip)+VPOCHA(L,J)
  195.             IEXI(K,IP)=1
  196.             III=III+1
  197.    9      CONTINUE
  198.    8    CONTINUE
  199.         SEGSUP,IPASS
  200.    5  CONTINUE
  201. *
  202.       IF (IIMPI.EQ.528) then
  203.         WRITE(IOIMP,1556)III,(DVAL(1,i),i=1,NPO)
  204.         WRITE(IOIMP,101)
  205.  1556   FORMAT(' III  DVAL ',I6,/,(1X,10E12.5))
  206.   101   FORMAT(' ON VIENT DE PASSER LA BOUCLE 5')
  207.       endif
  208. C
  209. C  **** ON REMPLIT DIN PAR LES VALEURS DE MCHPO1 POUR LES LAMBDAS
  210. C    mchpo1 est le vecteur initial FLX ( deplacement initial)
  211. C    din (icpr(i)) est le deplacement (FLX) initial du point i
  212. C
  213.       MCHPOI=MCHPO1
  214.       SEGACT,MCHPOI
  215.       DO 10 i=1,IPCHP(/1)
  216.         MSOUPO=IPCHP(I)
  217.         SEGACT,MSOUPO
  218.         MELEME=IGEOC
  219.         SEGACT,MELEME
  220.         MPOVAL=IPOVAL
  221.         SEGACT,MPOVAL
  222.         DO 11 J=1,NOCOMP(/2)
  223.           IF (NOCOMP(J).NE.'FLX ') GO TO 11
  224.           DO 12 K=1,NUM(/2)
  225.             IP=ICPR(NUM(1,K))
  226.             IF (IP.EQ.0) GO TO 12
  227.             DIN(IP)=VPOCHA(K,J)
  228.    12     CONTINUE
  229.    11   CONTINUE
  230.    10 CONTINUE
  231. *
  232.       IF (IIMPI.EQ.528) then
  233.         WRITE(IOIMP,102)
  234.         WRITE(IOIMP,666) (DIN(i),i=1,NPO)
  235.   102   FORMAT(' ON VIENT DE PASSER LA BOUCLE 10')
  236.   666   FORMAT('   DIN ',/,(1X,10E12.5))
  237.       ENDIF
  238. C
  239. C **** ON BOUCLE SUR LES RIGIDITE ET ON TESTE:
  240. C **** SI LE MULTI EXISTE DANS DVAL SON SIGNE PAR RAPPORT A IRIGEL(6,I)
  241. C **** SINON ON TESTE LA RELATION ELLE MEME A L'AIDE DES COEFF
  242. C **** DE LA MATRICE ET DE LA VALEUR DU LAMBDA INI (DANS DIN)
  243. C **** ON CREE UN TABLEAU CONTENANT LE NUMERO DES MATRICES A GARDER ET
  244. C **** LE NUMERO DE LA SOUS RIGIDITE
  245.       IPA=0
  246.       JG=NBO
  247.       SEGINI,MLENTI,trav1
  248.       DO 313 I=1,IRIGEL(/2)
  249.         ITY=IRIGEL(6,I)
  250.         IF (ITY.EQ.0) GO TO 313
  251.         MELEME=IRIGEL(1,I)
  252.         SEGACT,MELEME
  253.         DESCR=IRIGEL(3,I)
  254.         SEGACT,DESCR
  255.         DO 314 J=1,LISINC(/2)
  256.           IF(LISINC(J).EQ.'LX  ') GO TO 315
  257.   314   CONTINUE
  258.         CALL ERREUR(5)
  259.         RETURN
  260.   315   CONTINUE
  261.         xmatri=irigel(4,i)
  262.         segact,xmatri
  263.         JJ=NOELEP(J)
  264.         ITEMP=LISINC(/2)
  265.         SEGINI,IPASS
  266.         DO 316 J=1,ITEMP
  267.           DO 317 K=1,NIN
  268.             IF(LISINC(J).NE.NOMIN(K)) GO TO 317
  269.             IPASS(J)=K
  270.             GO TO 316
  271.   317     CONTINUE
  272.           CALL ERREUR (5)
  273.           RETURN
  274.   316   CONTINUE
  275. *  RECHERCHE DU MAX DES LAMDAS POUR LE CRITERE DE DECOLLEMENT D'APPUI
  276. *  RECHERCHE DU MAX DES DEPLACEMENT POUR LE CRITERE DE TRAVERSEE D'APPUI
  277. * ymcrit est le deplacement maxi de tous les points en relation unilateral
  278. * xmcrit est le maxi des multiplicateurs existant dans le chpoin de depla
  279. * remarque : ce chpoint de depla contient les multiplicateurs de contact
  280. * (pression)
  281.         DO 30 J=1,NUM(/2)
  282.          iel=j
  283.           IP=ICPR(NUM(JJ,J))
  284.           DO 31 K=1,ITEMP
  285.             IF(LISINC(K).EQ.'LX  ') GO TO 31
  286.             IPPP=NUM(NOELEP(K),J)
  287.             IPP=ICPR(IPPP)
  288. *  deplacement
  289.             YMCRIT=MAX(YMCRIT,ABS(DVAL(IPASS(K),IPP)))
  290.    31     CONTINUE
  291. *  jeu
  292.           YMCRIT=MAX(YMCRIT,ABS(DIN(IP)))
  293.           IF (IEXI(1,IP).EQ.0) GOTO 30
  294.           XMCRIT=MAX(XMCRIT,ABS(DVAL(1,IP)))
  295.   30    CONTINUE
  296. *       write (6,*) ' xmcrit ymcrit apres 30 ',xmcrit,ymcrit
  297.         SEGSUP,IPASS
  298. * on rajoute dans ymcrit la dimension de l'element
  299. **     write (6,*) ' avant 32 ymcrit num(/2) ',ymcrit,num(/2),num(/1),
  300. **   > re(/1),re(/3)
  301. ** uniquement si on travaille sur les deplacements
  302.         if (lisinc(2)(1:1).EQ.'U') then
  303.         idim1=idim+1
  304.         xcr=0.d0
  305.         do 32 iel=1,num(/2)
  306.           do 33 k=2,num(/1)-1
  307. *  le noeud 1 est le multiplicateur de lagrange
  308.              ip1=num(k,iel)
  309.              ip2=num(k+1,iel)
  310.              xp1=xcoor((ip1-1)*idim1+1)
  311.              yp1=xcoor((ip1-1)*idim1+2)
  312.              zp1=0.d0
  313.              if(idim.eq.3) zp1=xcoor((ip1-1)*idim1+3)
  314.              xp2=xcoor((ip2-1)*idim1+1)
  315.              yp2=xcoor((ip2-1)*idim1+2)
  316.              zp2=0.d0
  317.              if(idim.eq.3) zp2=xcoor((ip2-1)*idim1+3)
  318.              xcr=max(xcr,(xp2-xp1)**2+(yp2-yp1)**2+(zp2-zp1)**2)
  319.   33      continue
  320.   32    continue
  321.         ymcrit=max(ymcrit,sqrt(xcr))
  322.         endif
  323.  313  CONTINUE
  324. *
  325.       if (mchpo3.ne.0) then
  326. *        write (6,*) ' excit1 xmcrit avant ',xmcrit
  327.         do ip=1,ycpr(/1)
  328. ****      xmcrit=max(xmcrit,abs(ycpr(ip)))
  329.         enddo
  330. *        write (6,*) ' excit1 xmcrit apres ',xmcrit
  331.       endif
  332.  
  333.       XMCRIT=1D-9  *XMCRIT
  334.       YMCRIT=1D-9  *YMCRIT
  335. *     write (6,*) ' xmcrit ymcrit ',xmcrit,ymcrit
  336. * Critere en cas de frottement
  337.       yfcrit = YMCRIT
  338.  
  339. * Strategie lente plus laxiste
  340.       if (ityp.eq.3) then
  341.        xmcrit = xmcrit * 1d3
  342.        ymcrit = ymcrit * 1d3
  343.        yfcrit = yfcrit * 1d3
  344.       endif
  345. *
  346.       DO 13 I=1,IRIGEL(/2)
  347.         coer=coerig(i)
  348.         ITY=IRIGEL(6,I)
  349.         IF (ITY.EQ.0) GO TO 13
  350.         MELEME=IRIGEL(1,I)
  351.         SEGACT,MELEME
  352.         DESCR=IRIGEL(3,I)
  353.         SEGACT,DESCR
  354.         xMATRI=IRIGEL(4,I)
  355.         SEGACT,xMATRI
  356.         ITEMP=LISINC(/2)
  357.         SEGINI,IPASS
  358.         DO 14 J=1,ITEMP
  359.           IF (LISINC(J).EQ.'LX  ') GO TO 15
  360.    14   CONTINUE
  361.         CALL ERREUR(5)
  362.         RETURN
  363.    15   CONTINUE
  364.         JJ=NOELEP(J)
  365.         DO 16 J=1,ITEMP
  366.           DO 17 K=1,NIN
  367.             IF(LISINC(J).NE.NOMIN(K)) GO TO 17
  368.             IPASS(J)=K
  369.             GO TO 16
  370.    17     CONTINUE
  371.           CALL ERREUR (5)
  372.           RETURN
  373.    16   CONTINUE
  374.         DO 18 J=1,NUM(/2)
  375.           IPA=IPA+1
  376.           IPT=NUM(JJ,J)
  377.           IP=ICPR(IPT)
  378.           ityz=0
  379. *  eliminer condition redondantes pour le pas d'apres
  380.           iclred = .false. 
  381.           if (kcpr.ne.0) then
  382.             if (kcpr(num(noelep(1),j)).ne.0) then
  383. *  la condition redondante ayant ete augmentee on pourrait faire comme si elle n'etait pas presente
  384. *  mais en fait les criteres en reaction et en deplacement sont equivalents dans ce cas
  385.               if(iexi(1,ip).ne.0) iav(ipa)=2
  386.               iclred = .true. 
  387. ********      iexi(1,ip)=0
  388. *      write(6,*) 'condition redondante',ipa
  389.             endif
  390.           endif
  391. *
  392.             SS=0.D0
  393.             remax=0.d0
  394.             r_z=0.d0
  395.             DO 19 K=2,ITEMP
  396. ***           IF (LISINC(K).EQ.'LX  ') GO TO 19
  397.               if (ipass(k).eq.0) goto 19
  398.               IPPP=NUM(NOELEP(K),J)
  399.               IPP=ICPR(IPPP)
  400.               if (ipp.eq.0) goto 19
  401. *             write (6,*) ' k dval re ss',dval(ipass(k),ipp),
  402. *    >           re(1,k,j),ss
  403.               remax=max(remax,abs(re(1,k,j)*coer))
  404.               Sinc = DVAL(IPASS(K),IPP) * RE(1,k,j) * coer
  405.               r_z=max(r_z,abs(sinc))
  406.               SS = Sinc + SS
  407.    19       CONTINUE
  408.             IF (IIMPI.EQ.528) WRITE(IOIMP,529) IPPP,SS
  409. 529   FORMAT( ' LIBRE  ',I5,2X,E12.5)
  410.             r_z  = max(ABS(din(ip)),r_z)*1d-10
  411. *           write (6,*) 'r_z ymcrit ss',r_z,ymcrit,remax,ss
  412.             r_p1 = DIN(IP) + r_z
  413.             r_m1 = DIN(IP) - r_z
  414.             viol(ipa)=ss-din(ip)
  415. C
  416. C  ** CAS OU LE BLOQUAGE N'ETAIT PAS SOLLICITE
  417. C
  418.           IF (IEXI(1,IP).EQ.0.or.iclred) THEN 
  419. CC
  420.  
  421. *  cas du frottement
  422.             if (ity.eq.2) then
  423.               ityz=jcpr(ipt)
  424. *              write(ioimp,*) 'ipt jcpr ycpr ',ipt,jcpr(ipt),ycpr(ipt)
  425.               if (ityz.ne.0) ityz=sign(1.1D0,ycpr(ipt))
  426. *              if (ityz.eq.0) write (6,*) ' frottement -1 tyz ',ityz
  427. *              write(ioimp,*) '1 dans excite ityz ',ityz
  428. *  apparamment il faut etre plus laxiste pour le frottement
  429. *  peut etre pas finalement
  430.               if (ityz.eq.0) goto 20
  431.               IF (ITYz.EQ. 1.AND.SS.LE.r_p1+yfcrit*remax) GOTO 20
  432.               IF (ITYz.EQ.-1.AND.SS.GE.r_m1-yfcrit*remax) GOTO 20
  433.             endif
  434.             IF (ITY.EQ. 1.AND.SS.LE.r_p1+YMCRIT*remax) GOTO 20
  435.             IF (ITY.EQ.-1.AND.SS.GE.r_m1-YMCRIT*remax) GOTO 20
  436.             LECT(IPA)=ITY
  437.             IF(ity.eq.2) lect(ipa)=ityz*2
  438.        if (iimpi.eq.1967)
  439.      >     write (6,*) ' ss din ymcrit nouveau blocage '
  440.      $   ,ss,din(ip),ymcrit,ipa,viol(ipa),ipt,ityz,ityp
  441. *     on a un (1) nouveau blocage on arrete
  442.    20       CONTINUE
  443. C
  444. C  ** CAS OU LE BLOQUAGE ETAIT SOLLICITE  PETIT PROBLEME DE TEST DE
  445. C     PRECISION SUR LA VALEUR DE LA REACTION
  446. C
  447.           ELSE
  448.             iav(ipa)=1
  449.             SS=DVAL(1,IP)
  450.             viol(ipa)=ss
  451. *           write (6,*) ' ss xmcrit ',ss,xmcrit
  452.             IF(IIMPI.EQ.528) WRITE(IOIMP,530) NUM(JJ,J),SS
  453.   530 FORMAT('  BLOQUER ' ,I5,2X,E12.5)
  454. *  cas du frottement
  455.             if (ity.eq.2) then
  456.               ityz=jcpr(ipt)
  457.               if (ityz.ne.0) ityz=sign(1.1D0,ycpr(ipt))
  458. *              if (ityz.eq.0) write (6,*) ' frottement -2 ityz ',ityz
  459. *  apparamment il faut etre plus laxiste pour le frottement
  460. *  peut etre pas finalement
  461.               if (ityz.eq.0) goto 21
  462.               IF(ITYz.EQ.1.AND.SS.LT. -XMCRIT) GO TO 21
  463.               IF(ITYz.EQ.-1.AND.SS.GT.+XMCRIT) GO TO 21
  464.             endif
  465.             IF(ITY.EQ.1.AND.SS.LT. -XMCRIT) GO TO 21
  466.             IF(ITY.EQ.-1.AND.SS.GT.+XMCRIT) GO TO 21
  467.             LECT(IPA)=ITY
  468.             IF(ity.eq.2) lect(ipa)=ityz*2
  469.             goto 18
  470.    21       CONTINUE
  471. C      write (6,*) ' ss xmcrit ',ss,xmcrit,ityz,i,ipa,
  472. C     >  num(3,j),num(4,j),num(5,j)
  473.       if (iimpi.eq.1967)
  474.      > write (6,*) ' appui disparait '
  475.      $   ,ss,din(ip),ymcrit,ipa,viol(ipa),ipt,ityz,ityp
  476.           ENDIF
  477.    18   CONTINUE
  478.         SEGSUP,IPASS
  479.         SEGDES,DESCR,xMATRI
  480.    13 CONTINUE
  481.       IF(IIMPI.EQ.528) WRITE(IOIMP,*) 'ON VIENT DE PASSER LA BOUCLE 13'
  482. C
  483. C  **** DANS LECT ON DIT SI LA JEEME RIGI ELEMTAIRE EST A CONSERVER
  484. C
  485. *  on prevoit la place pour mettre a la fin les conditions redondantes
  486.       NRIGEL=IRIGEL(/2)*2
  487.       SEGINI,RI2
  488.       RI2.MTYMAT=mrigid.MTYMAT
  489.       RI2.IFORIG=mrigid.IFORIG
  490. C
  491. C  **** CAS OU IL N'Y A RIEN A GARDER  ON CREE UNE RIGIDITE VIDE
  492. C  **** POUR POUVOIR CONTINUER D'ITERER SI IL Y A LIEU
  493. C
  494.       IF (NRIGEL.EQ.0) THEN
  495.         IF (IIMPI.EQ.528) WRITE(IOIMP,*) ' IL N''Y A RIEN A CREER'
  496.         GO TO 50
  497.       ENDIF
  498. *
  499. *** ne garder que ce qui viole de plus de xcrot du max
  500. *    desactive pour le moment
  501. *
  502. * recherche du max
  503.       violmf=-1.
  504.       violmd=-1.
  505.       imf=0
  506.       imd=0
  507.       do 40 i=1,nbo
  508.         violm=abs(viol(i))
  509.         if (iav(i).eq.0.and.lect(i).ne.0) then
  510.           if (violm.gt.violmd) then
  511.             violmd=violm
  512.             imd=i
  513.           endif
  514.         else if (iav(i).ne.0.and.lect(i).eq.0) then
  515.           if (violm.gt.violmf) then
  516.             violmf=violm
  517.             imf=i
  518.           endif
  519.         endif
  520.   40  continue
  521.         rvd = 0.30*violmd
  522.         rvf = 0.90*violmf
  523.         if (ityp.eq.3) then
  524.         rvd = 0.9999*violmd
  525.         rvf = 0.9999*violmf
  526.         endif
  527.         do 41 i=1,nbo
  528.           violm=abs(viol(i))
  529.            if (iav(i).eq.0.and.lect(i).ne.0.and.violm.lt.rvd)
  530.      >      lect(i)=0
  531.   41    continue
  532. *  passe d'introduction de bruit en cas de conditions redondantes:
  533. *  on en elimine une en changeant laquelle a chaque appel
  534. ***     if(ityp.eq.3) then
  535.         if(kcpr.ne.0) then
  536.         icntr=mod(icntr,nbo)+1
  537.         icnt=0
  538.         do i=1,nbo
  539.           if (lect(i).ne.0) then
  540.             icnt=icnt+1
  541.               if (icnt.eq.icntr) then
  542.                 lect(i)=sign(4,lect(i))
  543.                 iav(i)=1
  544. ***               write(6,*) ' condition relachee  pour ',i
  545.                 goto 60
  546.               endif
  547.           endif
  548.         enddo
  549.   60    continue
  550.         icnts=mod(icnts,nbo)+1
  551.         icns=0
  552.         do i=1,nbo
  553.           if (lect(i).ne.0.and.iav(i).eq.2) then
  554.             icns=icns+1
  555.               if (icns.eq.icnts) then
  556.                 lect(i)=sign(4,lect(i))
  557.                 iav(i)=1
  558. ***               write(6,*) ' condition redondante relachee  pour ',i
  559.                 goto 62
  560.               endif
  561.           endif
  562.         enddo
  563.   62    continue
  564.         endif
  565.  
  566. C
  567. C  ** IL FAUT CREER UNE RIGIDITE
  568. C
  569.       IPA= 0
  570.       IRI=0
  571. *  premiere passe pour mettre les conditions redondantes
  572. *  en tete
  573.       DO 25 I =1,IRIGEL(/2)
  574.         IF(IRIGEL(6,I).EQ.0) GO TO 25
  575.         MELEME=IRIGEL(1,I)
  576.         SEGACT,MELEME
  577.         nbelei=num(/2)
  578.         NELRIG=0
  579.         DO 27 J=1,nbelei
  580.           IF(LECT(IPA+J).EQ.0) goto 27
  581.           IF(abs(LECT(IPA+J)).EQ.4) goto 27
  582.            if (kcpr.eq.0) goto 27
  583.            if (kcpr(num(1,j)).eq.0) goto 27
  584.            NELRIG=NELRIG+1
  585.   27    CONTINUE
  586.         if (nelrig.eq.0) goto 26
  587. **      write(6,*) ' nelrig redondant ',nelrig
  588.         IRI=IRI+1
  589.         xMATRI=IRIGEL(4,I)
  590.         SEGACT,xMATRI
  591.         nligrp=re(/2)
  592.         nligrd=re(/1)
  593.         SEGINI,xMATR1
  594.         RI2.IRIGEL(4,IRI)=xMATR1
  595.         RI2.IRIGEL(3,IRI)=IRIGEL(3,I)
  596.         RI2.IRIGEL(5,IRI)=IRIGEL(5,I)
  597.         RI2.IRIGEL(6,IRI)=IRIGEL(6,I)
  598.         RI2.IRIGEL(7,IRI)=IRIGEL(7,I)
  599.         xmatr1.symre=ri2.irigel(7,iri)
  600.         RI2.COERIG(IRI)=COERIG(I)
  601.         NBNN  =NUM(/1)
  602.         NBELEM=NELRIG
  603.         NBSOUS=0
  604.         NBREF =0
  605.         SEGINI,IPT1
  606.         IPT1.ITYPEL=ITYPEL
  607.         RI2.IRIGEL(1,IRI)=IPT1
  608.         I2=0
  609.         DO 28 J=1,nbelei
  610.           IF(LECT(IPA+J).EQ.0) goto 28
  611.           IF(abs(LECT(IPA+J)).EQ.4) goto 28
  612.            if (kcpr.eq.0) goto 28
  613.            if (kcpr(num(1,j)).eq.0) goto 28
  614.           I2=I2+1
  615.           DO 29 K=1,NUM(/1)
  616.             IPT1.NUM(K,I2)=NUM(K,J)
  617.  29       CONTINUE
  618.           do io=1,nligrp
  619.             do iu=1,nligrd
  620.               xmatr1.re(iu,io,i2)=re(iu,io,j)
  621.             enddo
  622.           enddo
  623. **            write (6,*) ' excit1 augmentation ',ipa+j
  624. **            xmatr1.re(1,1,i2)=-0.5
  625. **        endif
  626.  
  627.  28     CONTINUE
  628.         SEGDES,xMATR1
  629.  26     CONTINUE
  630.         IPA=IPA+nbelei
  631.    25 CONTINUE
  632. *  deuxieme passe pour mettre les conditions non redondantes
  633. *  en tete
  634.       nrige = irigel(/2)
  635.       ipa=0
  636.       do 125 ir=nrige+1,2*nrige
  637.         i=ir-nrige
  638.         IF(IRIGEL(6,I).EQ.0) GO TO 125
  639.         MELEME=IRIGEL(1,I)
  640.         SEGACT,MELEME
  641.         nbelei=num(/2)
  642.         NELRIG=0
  643.         DO 127 J=1,nbelei
  644.           IF(LECT(IPA+J).EQ.0) goto 127
  645.           IF(abs(LECT(IPA+J)).EQ.4) goto 127
  646.            if (kcpr.ne.0) then
  647.              if (kcpr(num(1,j)).ne.0) then
  648.              goto 127
  649.            endif
  650.            endif
  651.            NELRIG=NELRIG+1
  652.   127   CONTINUE
  653. **      if(nelrig.ne.0) write(6,*)'excit1 ',nelrig,'matrices deplacees'
  654.         if (nelrig.eq.0) goto 126
  655.         IRI=IRI+1
  656.         xMATRI=IRIGEL(4,I)
  657.         SEGACT,xMATRI
  658.         nligrp=re(/2)
  659.         nligrd=re(/1)
  660.         SEGINI,xMATR1
  661.         RI2.IRIGEL(4,IRI)=xMATR1
  662.         RI2.IRIGEL(3,IRI)=IRIGEL(3,I)
  663.         RI2.IRIGEL(5,IRI)=IRIGEL(5,I)
  664.         RI2.IRIGEL(6,IRI)=IRIGEL(6,I)
  665.         RI2.IRIGEL(7,IRI)=IRIGEL(7,I)
  666.         xmatr1.symre=ri2.irigel(7,iri)
  667.         RI2.COERIG(IRI)=COERIG(I)
  668.         NBNN  =NUM(/1)
  669.         NBELEM=NELRIG
  670.         NBSOUS=0
  671.         NBREF =0
  672.         SEGINI,IPT1
  673.         IPT1.ITYPEL=ITYPEL
  674.         RI2.IRIGEL(1,IRI)=IPT1
  675.         I2=0
  676.         DO 128 J=1,nbelei
  677.           IF(LECT(IPA+J).EQ.0) GO TO 128
  678.           IF(abs(LECT(IPA+J)).EQ.4) GO TO 128
  679.            if (kcpr.ne.0) then
  680.              if (kcpr(num(1,j)).ne.0) then
  681.              goto 128
  682.            endif
  683.            endif
  684.           I2=I2+1
  685.           DO K=1,NUM(/1)
  686.             IPT1.NUM(K,I2)=NUM(K,J)
  687.           ENDDO
  688.           do io=1,nligrp
  689.             do iu=1,nligrd
  690.               xmatr1.re(iu,io,i2)=re(iu,io,j)
  691.             enddo
  692.           enddo
  693.  128    CONTINUE
  694.         SEGDES,xMATR1,xMATRI
  695.  126    CONTINUE
  696.         IPA=IPA+nbelei
  697.  125  CONTINUE
  698. *
  699.       if (iri.ne.ri2.irigel(/2)) then
  700.          nrigel=iri
  701.          segadj,ri2
  702.       endif
  703.  
  704.    50 CONTINUE
  705.       SEGDES,RI2,MRIGID
  706.  
  707. *  indice de retour
  708.       indic = .true.
  709.       do 55 i = 1, nbo
  710.        if (iav(i).ne.lect(i)) indic = .false.
  711.   55  continue
  712. *     do il=1,lect(/1),5
  713. *      write (6,*)'  mlenti ',(lect(ill),ill=il,min(lect(/1),il+4))
  714. *     enddo
  715. **    do il=1,lect(/1),5
  716. **     write (6,*)'  viol   ',(viol(ill),ill=il,min(lect(/1),il+4))
  717. **    enddo
  718.  
  719.       SEGSUP,SNOMIN,ICPR,ITRAV,trav1
  720.       if (mchpo3.ne.0) segsup,jcpr,ycpr
  721.       if (kcpr.ne.0) segsup,kcpr
  722.  
  723.       END
  724.  
  725.  
  726.  
  727.  
  728.  

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