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

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