C EXCIT1 SOURCE PV090527 26/04/30 21:15:30 12529 SUBROUTINE EXCIT1(MRIGID,MCHPO1,MCHPO2,MCHPO3,MLENTI,indic,RI2, & mchpo4,ITYP) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) * modif : declarer inactifs les blocages redondants s'ils sont dans * mchpo4. Si celui-ci est non vide, on enleve les blocages inclus * dedans et on s'en tient la. -INC PPARAM -INC CCOPTIO -INC CCREEL -INC CCGEOME -INC SMRIGID -INC SMCHPOI -INC SMELEME -INC SMCOORD -INC SMLENTI save icntr,icnts DATA ICNTR/0/ DATA ICNTS/0/ SEGMENT SNOMIN CHARACTER*(LOCOMP) NOMIN(0) ENDSEGMENT SEGMENT ISOU(IRIGEL(/2)) SEGMENT ICPR(nbpts) SEGMENT JCPR(nbpts) SEGMENT KCPR(nbpts) SEGMENT YCPR(nbpts) SEGMENT ITRAV real*8 DVAL(NIN,NPO),DIN(NPO) integer IEXI(NIN,NPO) endsegment SEGMENT IPASS(ITEMP) segment trav1 integer iav(nbo) real*8 viol(nbo) endsegment segment trav2 integer lcpr(npo) real*8 zcpr(npo) endsegment data idec/0/ C C *** ICPR REFERENCIE DVAL C *** NOMIN CONTIENT LES INCONNUES REFERENCES PAR LES RELATIONS C *** DVAL CONTIENT LE RESULTAT DES DEPLACEMENTS ET DES LAMBDA C *** DIN CONTIENT LES VALEURS INITIALES DU SECOND MEMBRE C logical indic,iclred C idimp1 = idim+1 segact,mcoord C XMCRIT = xpetit/xzprec YMCRIT = xpetit/xzprec icpr=0 jcpr=0 ycpr=0 kcpr=0 itrav=0 trav2=0 C CHPOINT pour traiter les conditions redondantes. C N'identifier que les noeuds associes aux valeurs les plus C significatives * write (6,*) 'excit1 mchpo4 ',mchpo4 if (mchpo4.ne.0) then segact,mchpo4 do isous=1,mchpo4.ipchp(/1) msoupo=mchpo4.ipchp(isous) segact msoupo * write(6,*) 'nocomp(/1) nocomp(1)',nocomp(/1),nocomp(1) if (nocomp(/2).eq.1) then if (nocomp(1).eq.'LX') then ipt8=igeoc segact ipt8 if (ipt8.itypel.ne.1) then call erreur(16) return endif if (ipt8.num(/2).ne.0) then if (kcpr.eq.0) segini,kcpr mpoval=ipoval segact mpoval C C Valeur maximale xmax=0.D0 do i=1,ipt8.num(/2) xmax=max(xmax,abs(vpocha(i,1))) enddo C C Identifier les noeuds d'interet xma2=0.1*xmax do i=1,ipt8.num(/2) ** if(abs(vpocha(i,1)).gt.xpetit/xzprec) then if(abs(vpocha(i,1)).gt.xma2) then kcpr(ipt8.num(1,i))=1 endif enddo endif endif endif enddo endif * * mchpo3 est le champ des LX sorti de excfro jcpr donne existence * et ycpr donne valeur * write (ioimp,*) ' mchpo3 dans excit1 ',mchpo3 if (mchpo3.ne.0) then segini,jcpr,ycpr segact,mchpo3 do 100 isoupo=1,mchpo3.ipchp(/1) msoupo=mchpo3.ipchp(isoupo) segact,msoupo do 110 ic=1,nocomp(/2) if (nocomp(ic).eq.'LX ') goto 120 110 continue goto 140 120 continue ipt2=igeoc mpoval=ipoval segact,ipt2,mpoval do 130 iel=1,ipt2.num(/2) impt=ipt2.num(1,iel) jcpr(impt)=1 ycpr(impt)=ycpr(impt)+vpocha(iel,ic) 130 continue 140 continue 100 continue endif C C** INITIALISATION DE NOMIN ET ICPR icpr repere les points appartenant C** a la matrice de blocage mini maxi ou ? C SEGINI,SNOMIN,ICPR NOMIN(**)='LX ' NPO=0 NBO=0 * SEGACT,MRIGID DO 1 I=1,IRIGEL(/2) ITY=IRIGEL(6,I) IF (ITY.EQ.0) GO TO 1 * cas du frottement : petite verification if (ity.eq.2 .and. mchpo3.eq.0) then call erreur(721) return endif MELEME=IRIGEL(1,I) SEGACT,MELEME NBO=NBO+NUM(/2) IF(IIMPI.EQ.137) WRITE(IOIMP,3765) NBO 3765 FORMAT(' NBO ',I5) DO J=1,NUM(/2) DO K=1,NUM(/1) impt=NUM(K,J) IF (ICPR(impt).EQ.0) THEN NPO=NPO+1 ICPR(impt)=NPO ENDIF enddo enddo DESCR=IRIGEL(3,I) SEGACT,DESCR DO 3 J=1,LISINC(/2) DO 4 K=1,NOMIN(/2) IF (NOMIN(K).EQ.LISINC(J)) GO TO 3 4 CONTINUE NOMIN(**)=LISINC(J) 3 CONTINUE SEGDES,DESCR 1 CONTINUE NIN=NOMIN(/2) C if(kcpr.ne.0) segini,trav2 C IF(IIMPI.EQ.528) THEN WRITE(ioimp,90) NPO WRITE(IOIMP,9876)( NOMIN(i),i=1,NIN) 90 FORMAT(' ON VIENT DE PASSER LA BOUCLE 1 NPO ',I5) 9876 FORMAT(' NOMIN ' ,10(A8,1X)) ENDIF C C **** ON REMPLIT DVAL AVEC LES VALEURS DE MCHPO2 C mchpo2 est le champ de deplacement propose par RESOU C dval(j,icpr(i)) est la j eme composante du deplacement du point i C et iexi(j,icpr(i))=1 C SEGINI,ITRAV c MCHPOI=MCHPO2 SEGACT,MCHPOI DO 5 I=1,IPCHP(/1) MSOUPO=IPCHP(I) SEGACT,MSOUPO MELEME=IGEOC SEGACT,MELEME ITEMP=NOCOMP(/2) SEGINI,IPASS DO 6 K=1,ITEMP DO 7 J=1,NIN IF (NOMIN(J).EQ.NOCOMP(K)) THEN IPASS(K)=J GO TO 6 ENDIF 7 CONTINUE 6 CONTINUE IF (IIMPI.EQ.528) WRITE(IOIMP,1555) (IPASS(KHU),KHU=1,ITEMP) 1555 FORMAT(' IPASS ' ,9I10) MPOVAL=IPOVAL SEGACT,MPOVAL III=0 DO 8 J=1,ITEMP K=IPASS(J) IF (K.EQ.0) GO TO 8 DO 9 L=1,NUM(/2) IP=ICPR(NUM(1,L)) IF (IP.EQ.0) GO TO 9 * if (k.eq.1.and.vpocha(l,j).eq.0.d0 ) then * write (6,*) 'LX nul dans excit1' * goto 9 * endif DVAL(K,IP)=dval(k,ip)+VPOCHA(L,J) IEXI(K,IP)=1 III=III+1 9 CONTINUE 8 CONTINUE SEGSUP,IPASS 5 CONTINUE * IF (IIMPI.EQ.528) then WRITE(IOIMP,1556)III,(DVAL(1,i),i=1,NPO) WRITE(IOIMP,101) 1556 FORMAT(' III DVAL ',I6,/,(1X,10E12.5)) 101 FORMAT(' ON VIENT DE PASSER LA BOUCLE 5') endif C C **** ON REMPLIT DIN PAR LES VALEURS DE MCHPO1 POUR LES LAMBDAS C mchpo1 est le vecteur initial FLX ( deplacement initial) C din (icpr(i)) est le deplacement (FLX) initial du point i C MCHPOI=MCHPO1 SEGACT,MCHPOI DO 10 i=1,IPCHP(/1) MSOUPO=IPCHP(I) SEGACT,MSOUPO MELEME=IGEOC SEGACT,MELEME MPOVAL=IPOVAL SEGACT,MPOVAL DO 11 J=1,NOCOMP(/2) IF (NOCOMP(J).NE.'FLX ') GO TO 11 DO 12 K=1,NUM(/2) IP=ICPR(NUM(1,K)) IF (IP.EQ.0) GO TO 12 DIN(IP)=VPOCHA(K,J) 12 CONTINUE 11 CONTINUE 10 CONTINUE * IF (IIMPI.EQ.528) then WRITE(IOIMP,102) WRITE(IOIMP,666) (DIN(i),i=1,NPO) 102 FORMAT(' ON VIENT DE PASSER LA BOUCLE 10') 666 FORMAT(' DIN ',/,(1X,10E12.5)) ENDIF C C **** ON BOUCLE SUR LES RIGIDITE ET ON TESTE: C **** SI LE MULTI EXISTE DANS DVAL SON SIGNE PAR RAPPORT A IRIGEL(6,I) C **** SINON ON TESTE LA RELATION ELLE MEME A L'AIDE DES COEFF C **** DE LA MATRICE ET DE LA VALEUR DU LAMBDA INI (DANS DIN) C **** ON CREE UN TABLEAU CONTENANT LE NUMERO DES MATRICES A GARDER ET C **** LE NUMERO DE LA SOUS RIGIDITE IPA=0 JG=NBO SEGINI,MLENTI,trav1 DO 313 I=1,IRIGEL(/2) ITY=IRIGEL(6,I) IF (ITY.EQ.0) GO TO 313 MELEME=IRIGEL(1,I) SEGACT,MELEME DESCR=IRIGEL(3,I) SEGACT,DESCR DO 314 J=1,LISINC(/2) IF(LISINC(J).EQ.'LX ') GO TO 315 314 CONTINUE CALL ERREUR(5) RETURN 315 CONTINUE xmatri=irigel(4,i) segact,xmatri JJ=NOELEP(J) ITEMP=LISINC(/2) SEGINI,IPASS DO 316 J=1,ITEMP DO 317 K=1,NIN IF(LISINC(J).NE.NOMIN(K)) GO TO 317 IPASS(J)=K GO TO 316 317 CONTINUE CALL ERREUR (5) RETURN 316 CONTINUE * RECHERCHE DU MAX DES LAMDAS POUR LE CRITERE DE DECOLLEMENT D'APPUI * RECHERCHE DU MAX DES DEPLACEMENT POUR LE CRITERE DE TRAVERSEE D'APPUI * ymcrit est le deplacement maxi de tous les points en relation unilateral * xmcrit est le maxi des multiplicateurs existant dans le chpoin de depla * remarque : ce chpoint de depla contient les multiplicateurs de contact * (pression) DO 30 J=1,NUM(/2) iel=j IP=ICPR(NUM(JJ,J)) DO 31 K=1,ITEMP IF(LISINC(K).EQ.'LX ') GO TO 31 IPPP=NUM(NOELEP(K),J) IPP=ICPR(IPPP) * deplacement YMCRIT=MAX(YMCRIT,ABS(DVAL(IPASS(K),IPP))) 31 CONTINUE * jeu YMCRIT=MAX(YMCRIT,ABS(DIN(IP))) IF (IEXI(1,IP).EQ.0) GOTO 30 XMCRIT=MAX(XMCRIT,ABS(DVAL(1,IP))) 30 CONTINUE * write (6,*) ' xmcrit ymcrit apres 30 ',xmcrit,ymcrit SEGSUP,IPASS * on rajoute dans ymcrit la dimension de l'element ** write (6,*) ' avant 32 ymcrit num(/2) ',ymcrit,num(/2),num(/1), ** > re(/1),re(/3) ** uniquement si on travaille sur les deplacements if (lisinc(2)(1:1).EQ.'U') then idim1=idim+1 xcr=0.d0 do 32 iel=1,num(/2) do 33 k=2,num(/1)-1 * le noeud 1 est le multiplicateur de lagrange ip1=num(k,iel) ip2=num(k+1,iel) xp1=xcoor((ip1-1)*idim1+1) yp1=xcoor((ip1-1)*idim1+2) zp1=0.d0 if(idim.eq.3) zp1=xcoor((ip1-1)*idim1+3) xp2=xcoor((ip2-1)*idim1+1) yp2=xcoor((ip2-1)*idim1+2) zp2=0.d0 if(idim.eq.3) zp2=xcoor((ip2-1)*idim1+3) xcr=max(xcr,(xp2-xp1)**2+(yp2-yp1)**2+(zp2-zp1)**2) 33 continue 32 continue ymcrit=max(ymcrit,sqrt(xcr)) endif 313 CONTINUE * if (mchpo3.ne.0) then * write (6,*) ' excit1 xmcrit avant ',xmcrit do ip=1,ycpr(/1) **** xmcrit=max(xmcrit,abs(ycpr(ip))) enddo * write (6,*) ' excit1 xmcrit apres ',xmcrit endif XMCRIT=1D-9 *XMCRIT YMCRIT=1D-9 *YMCRIT * write (6,*) ' xmcrit ymcrit ',xmcrit,ymcrit * Critere en cas de frottement yfcrit = YMCRIT * Strategie lente plus laxiste if (ityp.eq.3) then xmcrit = xmcrit * 1d3 ymcrit = ymcrit * 1d3 yfcrit = yfcrit * 1d3 endif * * Conditions redondantes xxxx=-1.D0 imuok=0 * DO 13 I=1,IRIGEL(/2) coer=coerig(i) ITY=IRIGEL(6,I) IF (ITY.EQ.0) GO TO 13 MELEME=IRIGEL(1,I) SEGACT,MELEME DESCR=IRIGEL(3,I) SEGACT,DESCR xMATRI=IRIGEL(4,I) SEGACT,xMATRI ITEMP=LISINC(/2) SEGINI,IPASS DO 14 J=1,ITEMP IF (LISINC(J).EQ.'LX ') GO TO 15 14 CONTINUE CALL ERREUR(5) RETURN 15 CONTINUE JJ=NOELEP(J) DO 16 J=1,ITEMP DO 17 K=1,NIN IF(LISINC(J).NE.NOMIN(K)) GO TO 17 IPASS(J)=K GO TO 16 17 CONTINUE CALL ERREUR (5) RETURN 16 CONTINUE DO 18 J=1,NUM(/2) IPA=IPA+1 IPT=NUM(JJ,J) IP=ICPR(IPT) ityz=0 C C Traitement particulier pour les conditions redondantes iclred = .false. if (kcpr.ne.0) then if (kcpr(num(noelep(1),j)).ne.0) then * la condition redondante ayant ete augmentee on pourrait faire comme si elle n'etait pas presente * mais en fait les criteres en reaction et en deplacement sont equivalents dans ce cas ** if(iexi(1,ip).ne.0) iav(ipa)=2 iclred = .true. ******** iexi(1,ip)=0 * write(6,*) 'condition redondante',ipa endif endif * SS=0.D0 remax=0.d0 r_z=0.d0 DO 19 K=2,ITEMP *** IF (LISINC(K).EQ.'LX ') GO TO 19 if (ipass(k).eq.0) goto 19 IPPP=NUM(NOELEP(K),J) IPP=ICPR(IPPP) if (ipp.eq.0) goto 19 * write (6,*) ' k dval re ss',dval(ipass(k),ipp), * > re(1,k,j),ss remax=max(remax,abs(re(1,k,j)*coer)) Sinc = DVAL(IPASS(K),IPP) * RE(1,k,j) * coer r_z=max(r_z,abs(sinc)) SS = Sinc + SS 19 CONTINUE IF (IIMPI.EQ.528) WRITE(IOIMP,529) IPPP,SS 529 FORMAT( ' LIBRE ',I5,2X,E12.5) r_z = max(ABS(din(ip)),r_z)*1d-10 * write (6,*) 'r_z ymcrit ss',r_z,ymcrit,remax,ss r_p1 = DIN(IP) + r_z r_m1 = DIN(IP) - r_z viol(ipa)=ss-din(ip) C C ** CAS OU LE BLOQUAGE N'ETAIT PAS SOLLICITE C IF (IEXI(1,IP).EQ.0.or.iclred) THEN C C Condition de frottement if (ity.eq.2) then ityz=jcpr(ipt) * write(ioimp,*) 'ipt jcpr ycpr ',ipt,jcpr(ipt),ycpr(ipt) if (ityz.ne.0) ityz=sign(1.1D0,ycpr(ipt)) * if (ityz.eq.0) write (6,*) ' frottement -1 tyz ',ityz * write(ioimp,*) '1 dans excite ityz ',ityz * apparamment il faut etre plus laxiste pour le frottement * peut etre pas finalement if (ityz.eq.0) goto 20 YCRIT=YFCRIT else ityz = ity YCRIT=YMCRIT endif IF (ITYz.EQ. 1.AND.SS.LE.r_p1+YCRIT*remax) GOTO 20 IF (ITYz.EQ.-1.AND.SS.GE.r_m1-YCRIT*remax) GOTO 20 C C Nouvelle condition a prendre en compte LECT(IPA)=ITY IF(ity.eq.2) lect(ipa)=ityz*2 * if (iimpi.eq.1967) > write (6,*) ' ss din ymcrit nouveau blocage ' $ ,ss,din(ip),ymcrit,ipa,viol(ipa),ipt,ityz,ityp * on a un (1) nouveau blocage on arrete C 20 CONTINUE * IF (iclred) then zcpr(ip)=viol(ipa) xxxx = max(xxxx,abs(viol(ipa))) ENDIF ENDIF C C ** CAS OU LE BLOQUAGE ETAIT SOLLICITE C PETIT PROBLEME DE TEST DE PRECISION SUR LA VALEUR DE LA REACTION C IF (IEXI(1,IP).NE.0.or.iclred) THEN iav(ipa)=1 SS=DVAL(1,IP) viol(ipa)=ss * write (6,*) ' ss xmcrit ',ss,xmcrit IF(IIMPI.EQ.528) WRITE(IOIMP,530) NUM(JJ,J),SS 530 FORMAT(' BLOQUER ' ,I5,2X,E12.5) C C Condition de frottement if (ity.eq.2) then ityz=jcpr(ipt) if (ityz.ne.0) ityz=sign(1.1D0,ycpr(ipt)) * if (ityz.eq.0) write (6,*) ' frottement -2 ityz ',ityz * apparamment il faut etre plus laxiste pour le frottement * peut etre pas finalement if (ityz.eq.0) goto 21 else ITYz=ITY endif IF(ITYz.EQ. 1.AND.SS.LT.-XMCRIT) GOTO 21 IF(ITYz.EQ.-1.AND.SS.GT.+XMCRIT) GOTO 21 C C Conserver la condition LECT(IPA)=ITY IF(ity.eq.2) lect(ipa)=ityz*2 C IF (iclred) then imuok=imuok+1 lcpr(ip)=1 ENDIF C GOTO 18 C 21 CONTINUE C write (6,*) ' ss xmcrit ',ss,xmcrit,ityz,i,ipa, C > num(3,j),num(4,j),num(5,j) if (iimpi.eq.1967) > write (6,*) ' appui disparait ' $ ,ss,din(ip),ymcrit,ipa,viol(ipa),ipt,ityz,ityp ENDIF C 18 CONTINUE SEGSUP,IPASS SEGDES,DESCR,xMATRI 13 CONTINUE IF(IIMPI.EQ.528) WRITE(IOIMP,*) 'ON VIENT DE PASSER LA BOUCLE 13' C C **** DANS LECT ON DIT SI LA JEEME RIGI ELEMTAIRE EST A CONSERVER C NRIGEL=IRIGEL(/2)*2 SEGINI,RI2 RI2.MTYMAT=mrigid.MTYMAT RI2.IFORIG=mrigid.IFORIG C C **** CAS OU IL N'Y A RIEN A GARDER ON CREE UNE RIGIDITE VIDE C **** POUR POUVOIR CONTINUER D'ITERER SI IL Y A LIEU C IF (NRIGEL.EQ.0) THEN IF (IIMPI.EQ.528) WRITE(IOIMP,*) ' IL N''Y A RIEN A CREER' GO TO 50 ENDIF * *** ne garder que ce qui viole de plus de xcrot du max * desactive pour le moment * * recherche du max violmf=-1. violmd=-1. imf=0 imd=0 do 40 i=1,nbo violm=abs(viol(i)) if (iav(i).eq.0.and.lect(i).ne.0) then if (violm.gt.violmd) then violmd=violm imd=i endif else if (iav(i).ne.0.and.lect(i).eq.0) then if (violm.gt.violmf) then violmf=violm imf=i endif endif 40 continue rvd = 0.30*violmd rvf = 0.90*violmf xxx = 0.90*xxxx if (ityp.eq.3) then rvd = 0.9999*violmd rvf = 0.9999*violmf xxx = 0.9999*xxxx endif do 41 i=1,nbo violm=abs(viol(i)) if (iav(i).eq.0.and.lect(i).ne.0.and.violm.lt.rvd) lect(i)=0 41 continue C C ** IL FAUT CREER UNE RIGIDITE C IPA=0 IRI=0 DO 25 I =1,IRIGEL(/2) IF(IRIGEL(6,I).EQ.0) GOTO 25 MELEME=IRIGEL(1,I) SEGACT,MELEME nbelei=num(/2) NELRIG=0 DO 27 J=1,nbelei C C Traitement de conditions redondantes if (kcpr.ne.0) then ipt=num(1,j) itt=icpr(ipt) if (imuok.eq.0) then if (kcpr(ipt).ne.0) then if (abs(zcpr(itt)).ge.xxx) then lect(ipa+j)=ity if (ity.eq.2) then ityz=jcpr(ipt) if (ityz.ne.0) ityz=sign(1.1D0,ycpr(ipt)) lect(ipa+j)=ityz*2 ENDIF endif endif else if (kcpr(ipt).ne.0) then if (abs(zcpr(itt)).lt.xxx.and.imuok.gt.1) then lect(ipa+j)=0 endif endif endif endif C IF(LECT(IPA+J).EQ.0) goto 27 C NELRIG=NELRIG+1 27 CONTINUE if (nelrig.eq.0) goto 26 IRI=IRI+1 CC if (iri.gt.nrigel) iri=1 xMATRI=IRIGEL(4,I) SEGACT,xMATRI nligrp=re(/2) nligrd=re(/1) rigrel=0 SEGINI,xMATR1 RI2.IRIGEL(4,IRI)=xMATR1 RI2.IRIGEL(3,IRI)=IRIGEL(3,I) RI2.IRIGEL(5,IRI)=IRIGEL(5,I) RI2.IRIGEL(6,IRI)=IRIGEL(6,I) RI2.IRIGEL(7,IRI)=IRIGEL(7,I) xmatr1.symre=ri2.irigel(7,iri) RI2.COERIG(IRI)=COERIG(I) NBNN =NUM(/1) NBELEM=NELRIG NBSOUS=0 NBREF =0 SEGINI,IPT1 IPT1.ITYPEL=ITYPEL RI2.IRIGEL(1,IRI)=IPT1 I2=0 DO 28 J=1,nbelei IF(LECT(IPA+J).EQ.0) goto 28 I2=I2+1 DO 29 K=1,NUM(/1) IPT1.NUM(K,I2)=NUM(K,J) 29 CONTINUE do io=1,nligrp do iu=1,nligrd xmatr1.re(iu,io,i2)=re(iu,io,j) enddo enddo ** write (6,*) ' excit1 augmentation ',ipa+j ** xmatr1.re(1,1,i2)=-0.5 ** endif 28 CONTINUE SEGDES,xMATR1 26 CONTINUE IPA=IPA+nbelei 25 CONTINUE * if (iri.ne.ri2.irigel(/2)) then nrigel=iri segadj,ri2 endif C 50 CONTINUE SEGDES,RI2,MRIGID * indice de retour indic = .true. do 55 i = 1, nbo if (iav(i).ne.lect(i)) indic = .false. 55 continue * do il=1,lect(/1),5 * write (6,*)' mlenti ',(lect(ill),ill=il,min(lect(/1),il+4)) * enddo ** do il=1,lect(/1),5 ** write (6,*)' viol ',(viol(ill),ill=il,min(lect(/1),il+4)) ** enddo C SEGSUP,SNOMIN,ICPR,ITRAV,trav1 if (mchpo3.ne.0) segsup,jcpr,ycpr if (kcpr.ne.0) segsup,kcpr,trav2 END