C IMPO32 SOURCE PV 13/12/13 21:15:04 7883 * impo bloc en 3D SUBROUTINE IMPO32 IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC CCOPTIO -INC CCREEL -INC CCGEOME -INC SMCOORD -INC SMELEME -INC SMRIGID -INC SMCHPOI * Contact en formulation forte : * Segments utiles dans les cas dits pathologiques segment mfopa1 integer lpoin(nfopa1,3) real*8 xpoin(nfopa1,3) endsegment segment mfopa2 integer lsegm(nfopa2,4) real*8 xsegm(nfopa2,3) endsegment * directions moyennes aux sommets segment mfopa3 real*8 xms(nbpts),yms(nbpts),zms(nbpts) endsegment * Contact en formulation faible : * Segment utile pour le calcul de l'intersection des 2 triangles * projetes dans le plan de contact intermediaire segment mfaible real*8 cPT1C(3,3),cPT2C(3,3), cPT1A(4,2),cPT2A(4,2) real*8 cPIn0(6,2), cPIn(6,2), test(6) real*8 vT1A(3,2), vT2A(3,2) real*8 SuT1A, SuT2A, SuPIn *DBG-F real*8 xGIn,yGIn, b1T1,b2T1,b3T1, b1T2,b2T2,b3T2 endsegment PARAMETER ( X1s3=0.333333333333333333333333333333333333333333D0 ) character*4 modepl(3),moforc(3) data modepl /'UX ','UY ','UZ '/ data moforc /'FX ','FY ','FZ '/ * idimp1 = IDIM + 1 * * Lecture du maillage support des conditions de contact * call lirobj('MAILLAGE',ipt1,1,iretou) if (ierr.ne.0) return * * Activation du maillage et petites verifications * segact ipt1 nbno1 = ipt1.num(/1) nbel1 = ipt1.num(/2) if (ipt1.lisous(/1).ne.0) call erreur(25) if (ipt1.itypel.ne.22) call erreur(16) if (nbno1.ne.5 .and. nbno1.ne.7) call erreur(16) if (ierr.ne.0) goto 900 * * Increment sur le nombre d'elements de contact retenus et utilises * dans le chpoint (mpoval et igeoc) et rigidite pour agrandir les * longueurs des segments adequats incnel = 2000 * * Creation de la raideur des conditions de contact * Remplissage de l'entete commun * nrigel = 1 segini,mrigid ichole = 0 imgeo1 = 0 imgeo2 = 0 isupeq = 0 iforig = ifomod coerig(1) = 1. mtymat='RIGIDITE' * * Creation du chpoint de depi * nat=1 nsoupo=1 segini mchpoi mtypoi='DEPIMP' mochde='engendré par impose' ifopoi=ifomod jattri(1)=2 * nc=1 segini msoupo ipchp(1)=msoupo nocomp(1)='FLX ' nocons(1)=' ' * nbnn =1 nbelem=incnel nbref =0 nbsous=0 segini ipt7 ipt7.itypel=1 igeoc=ipt7 * n=incnel segini mpoval ipoval=mpoval * * Nombre de noeuds dans le chpoint (en totalite) : ipt7 et mpoval nelch = 0 *======================================================================= * Formulation "forte" (standard) du contact : *======================================================================= * Element de contact 3D a 5 noeuds if (nbno1.eq.7) goto 500 * Calcul des directions moyennes nbpts = xcoor(/1)/(idim+1) segini mfopa3 do 820 iel=1,nbel1 * ip1=ipt1.num(2,iel) ip2=ipt1.num(3,iel) ip3=ipt1.num(4,iel) ipv1 = (ip1-1)*idimp1 xp1 = xcoor(ipv1+1) yp1 = xcoor(ipv1+2) zp1 = xcoor(ipv1+3) ipv2 = (ip2-1)*idimp1 xp2 = xcoor(ipv2+1) yp2 = xcoor(ipv2+2) zp2 = xcoor(ipv2+3) ipv3 = (ip3-1)*idimp1 xp3 = xcoor(ipv3+1) yp3 = xcoor(ipv3+2) zp3 = xcoor(ipv3+3) * * normale au plan (123) * x12 = xp2 - xp1 y12 = yp2 - yp1 z12 = zp2 - zp1 x23 = xp3 - xp2 y23 = yp3 - yp2 z23 = zp3 - zp2 xn = y12*z23 - z12*y23 yn = z12*x23 - x12*z23 zn = x12*y23 - y12*x23 sn = sqrt (xn*xn + yn*yn + zn*zn) xn = xn/sn yn = yn/sn zn = zn/sn xms(ip1)=xms(ip1)+xn yms(ip1)=yms(ip1)+yn zms(ip1)=zms(ip1)+zn xms(ip2)=xms(ip2)+xn yms(ip2)=yms(ip2)+yn zms(ip2)=zms(ip2)+zn xms(ip3)=xms(ip3)+xn yms(ip3)=yms(ip3)+yn zms(ip3)=zms(ip3)+zn 820 continue do 821 ip = 1,nbpts sn = xms(ip)*xms(ip)+yms(ip)*yms(ip)+zms(ip)*zms(ip) sn = max(sqrt(abs(sn)),1d-35) xms(ip)=xms(ip)/sn yms(ip)=yms(ip)/sn zms(ip)=zms(ip)/sn 821 continue * * Nombre de "cas pathologiques" de type noeud-noeud (type 1) nelpa1 = 0 * Nombre de "cas pathologiques" de type noeud-segment (type 2) nelpa2 = 0 * Critere de detection des cas pathologiques dcrit = 1.d-2 * Segments lies aux cas pathologiques (types 1 et 2) nfopa1 = incnel nfopa2 = incnel segini,mfopa1,mfopa2 * * Relation du type 3 : noeud-triangle *------------------------------------- * * creation du meleme associe a la relation * nbsous=0 nbref =0 nbnn =5 nbelem=incnel segini meleme itypel=22 irigel(1,nrigel)=meleme * * creation du descriptif commun a toutes les raideurs * nligrp = 13 nligrd = nligrp segini,descr lisinc(1)='LX ' lisdua(1)='FLX ' noelep(1)=1 noeled(1)=1 do 10 i=2, nligrp, 3 lisinc(i) =modepl(1) lisinc(i+1)=modepl(2) lisinc(i+2)=modepl(3) lisdua(i) =moforc(1) lisdua(i+1)=moforc(2) lisdua(i+2)=moforc(3) noelep(i) =(i+5)/3 noelep(i+1)=noelep(i) noelep(i+2)=noelep(i) noeled(i) =noelep(i) noeled(i+1)=noelep(i) noeled(i+2)=noelep(i) 10 continue segdes,descr irigel(3,nrigel) = descr * * creation du segment xmatri * nelrig = incnel segini,xmatri irigel(4,nrigel) = xmatri * * ce qu'on cree est unilateral * irigel(6,nrigel)=1 * * ce qu'on cree est symetrique * irigel(7,nrigel)=0 * * Nombre d'elements dans la rigidite du type 3 nelri3 = 0 * Boucle sur le maillage support du contact(frottement) * ip1o = 0 ip2o = 0 ip3o = 0 DO 20 iel=1,nbel1 * ip1=ipt1.num(2,iel) ip2=ipt1.num(3,iel) ip3=ipt1.num(4,iel) if (ip1.ne.ip1o.or.ip2.ne.ip2o.or.ip3.ne.ip3o) then ip1o=ip1 ip2o=ip2 ip3o=ip3 * Recuperation des coordonees des noeuds du triangle ipv = (ip1-1)*idimp1 xp1 = xcoor(ipv+1) yp1 = xcoor(ipv+2) zp1 = xcoor(ipv+3) ipv = (ip2-1)*idimp1 xp2 = xcoor(ipv+1) yp2 = xcoor(ipv+2) zp2 = xcoor(ipv+3) ipv = (ip3-1)*idimp1 xp3 = xcoor(ipv+1) yp3 = xcoor(ipv+2) zp3 = xcoor(ipv+3) * Centre de gravite du triangle xg = (xp1+xp2+xp3) * X1s3 yg = (yp1+yp2+yp3) * X1s3 zg = (zp1+zp2+zp3) * X1s3 endif jp =ipt1.num(5,iel) ipv = (jp-1)*idimp1 xp = xcoor(ipv+1) yp = xcoor(ipv+2) zp = xcoor(ipv+3) * Verification que pas relation du point sur lui meme if (jp.eq.ip1) goto 20 if (jp.eq.ip2) goto 20 if (jp.eq.ip3) goto 20 * * Verification si autre point dans la zone de selection * verif distance au centre de gravite d1 = (xg-xp1)**2 + (yg-yp1)**2 + (zg-zp1)**2 d2 = (xg-xp2)**2 + (yg-yp2)**2 + (zg-zp2)**2 d3 = (xg-xp3)**2 + (yg-yp3)**2 + (zg-zp3)**2 dist = max(d1,d2,d3) ** dp = ((xg-xp)**2 + (yg-yp)**2 + (zg-zp)**2) * 0.25 dp = ((xg-xp)**2 + (yg-yp)**2 + (zg-zp)**2) * 0.64 c*OF A voir -> dp = (xg-xp)**2 + (yg-yp)**2 + (zg-zp)**2 if (dp.gt.dist*2) goto 20 C*DBG write(ioimp,*) 'contact test distance ok',dp,d1,d2,d3 dist = SQRT(dist) * verif position par rapport aux cotes * cote 1 2 x12 = xp2 - xp1 y12 = yp2 - yp1 z12 = zp2 - zp1 dv = sqrt((x12*x12)+(y12*y12)+(z12*z12)) if (.not.(dv.lt.1d-10) .and. .not.(dv.gt.1d-10)) & write(ioimp,*) 'Prob 1.1 - impo32' if (abs(dv).le.xpetit) write(ioimp,*) 'Prob 1.2 - impo32' xv = x12 / dv yv = y12 / dv zv = z12 / dv scal = (xp-xp1)*xv + (yp-yp1)*yv + (zp-zp1)*zv xm = xp1 + scal*xv ym = yp1 + scal*yv zm = zp1 + scal*zv scal = (xg-xm)*xv + (yg-ym)*yv + (zg-zm)*zv xn = (xg-xm) - scal*xv yn = (yg-ym) - scal*yv zn = (zg-zm) - scal*zv dn = sqrt(xn*xn+yn*yn+zn*zn) dpm = sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2) scal = (xp-xm)*xn + (yp-ym)*yn + (zp-zm)*zn if (dpm.gt.1d-3*dv .and. scal.lt.-0.8d0*dn*dpm) goto 20 C*DBG write(ioimp,*) 'contact test position 1 ok', C*DBG & scal/(dn*dpm),scal,dn,dpm * cote 2 3 x23 = xp3 - xp2 y23 = yp3 - yp2 z23 = zp3 - zp2 dv = sqrt((x23*x23)+(y23*y23)+(z23*z23)) if (.not.(dv.lt.1d-10) .and. .not.(dv.gt.1d-10)) & write(ioimp,*) 'Prob 2.1 - impo32' if (abs(dv).le.xpetit) write(ioimp,*) 'Prob 2.2 - impo32' xv = x23 / dv yv = y23 / dv zv = z23 / dv scal = (xp-xp2)*xv + (yp-yp2)*yv + (zp-zp2)*zv xm = xp2 + scal*xv ym = yp2 + scal*yv zm = zp2 + scal*zv scal = (xg-xm)*xv + (yg-ym)*yv + (zg-zm)*zv xn = (xg-xm) - scal*xv yn = (yg-ym) - scal*yv zn = (zg-zm) - scal*zv dn = sqrt(xn*xn+yn*yn+zn*zn) dpm = sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2) scal = (xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn if (dpm.gt.1d-3*dv .and. scal.lt.-0.8d0*dn*dpm) goto 20 C*DBG write(ioimp,*) 'contact test position 2 ok', C*DBG & scal/(dn*dpm),scal,dn,dpm * cote 3 1 x31 = xp1 - xp3 y31 = yp1 - yp3 z31 = zp1 - zp3 dv = sqrt((x31*x31)+(y31*y31)+(z31*z31)) if (.not.(dv.lt.1d-10) .and. .not.(dv.gt.1d-10)) & write(ioimp,*) 'Prob 3.1 - impo32' if (abs(dv).le.xpetit) write(ioimp,*) 'Prob 3.2 - impo32' xv = x31 / dv yv = y31 / dv zv = z31 / dv scal = (xp-xp3)*xv + (yp-yp3)*yv + (zp-zp3)*zv xm = xp3 + scal*xv ym = yp3 + scal*yv zm = zp3 + scal*zv scal = (xg-xm)*xv + (yg-ym)*yv + (zg-zm)*zv xn = (xg-xm) - scal*xv yn = (yg-ym) - scal*yv zn = (zg-zm) - scal*zv dn = sqrt(xn*xn+yn*yn+zn*zn) dpm = sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2) scal = (xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn if (dpm.gt.1d-3*dv .and .scal.lt.-0.8d0*dn*dpm) goto 20 C*DBG write(ioimp,*) 'contact test position 3 ok', C*DBG & scal/(dn*dpm),scal,dn,dpm * * on a un point ou imposer la relation C*DBG write(ioimp,*) 'contact potentiel ' * Quelle est la relation ??? * * direction de la relation = normale au plan (123) * xn = y12*z23 - z12*y23 yn = z12*x23 - x12*z23 zn = x12*y23 - y12*x23 dn = sqrt(xn*xn+yn*yn+zn*zn) if (.not.(dn.lt.1d-10).and. .not.(dn.gt.1d-10)) & write(ioimp,*) 'Prob 4.1 - impo32' if (abs(dn).le.xpetit) write(ioimp,*) 'Prob 4.2 - impo32' xn = xn / dn yn = yn / dn zn = zn / dn if (abs(xn).lt.1d-12) xn = 0.d0 if (abs(yn).lt.1d-12) yn = 0.d0 if (abs(zn).lt.1d-12) zn = 0.d0 * * Distance (jeu) du point jp au plan de contact * Puis calcul de la projection sur le plan de contact * xjeu = (xp-xg)*xn + (yp-yg)*yn + (zp-zg)*zn * verif bon cote (a la tolerance pres) * write (6,*) ' xjeu ',xjeu,dist if (xjeu.lt.-0.8d0*dist) goto 20 * * Recherche de ses coordonnées barycentriques * qui sont les surfaces des sous triangles * xb1 = (yp-yp2)*(zp-zp3)-(zp-zp2)*(yp-yp3) yb1 = (zp-zp2)*(xp-xp3)-(xp-xp2)*(zp-zp3) zb1 = (xp-xp2)*(yp-yp3)-(yp-yp2)*(xp-xp3) xb2 = (yp-yp3)*(zp-zp1)-(zp-zp3)*(yp-yp1) yb2 = (zp-zp3)*(xp-xp1)-(xp-xp3)*(zp-zp1) zb2 = (xp-xp3)*(yp-yp1)-(yp-yp3)*(xp-xp1) xb3 = (yp-yp1)*(zp-zp2)-(zp-zp1)*(yp-yp2) yb3 = (zp-zp1)*(xp-xp2)-(xp-xp1)*(zp-zp2) zb3 = (xp-xp1)*(yp-yp2)-(yp-yp1)*(xp-xp2) b1 = xb1*xn + yb1*yn + zb1*zn b2 = xb2*xn + yb2*yn + zb2*zn b3 = xb3*xn + yb3*yn + zb3*zn bt = b1 + b2 + b3 if (.not.(bt.lt.1d-10) .and. .not.(bt.gt.1d-10)) & write(ioimp,*) 'Prob 5.1 - impo32' if (abs(bt).le.xpetit) write(ioimp,*) 'Prob 5.2 - impo32' b1 = b1 / bt b2 = b2 / bt b3 = b3 / bt b1c = b1 b2c = b2 b3c = b3 * recalcul de la direction et du jeu en fonction des normales aux sommets if (b3c.lt.0) then b3c = 0.d0 b1c = (xp-xp2)*(xp1-xp2)+(yp-yp2)*(yp1-yp2)+(zp-zp2)*(zp1-zp2) b2c = (xp-xp1)*(xp2-xp1)+(yp-yp1)*(yp2-yp1)+(zp-zp1)*(zp2-zp1) elseif (b1c.lt.0.d0) then b1c = 0.d0 b2c = (xp-xp3)*(xp2-xp3)+(yp-yp3)*(yp2-yp3)+(zp-zp3)*(zp2-zp3) b3c = (xp-xp2)*(xp3-xp2)+(yp-yp2)*(yp3-yp2)+(zp-zp2)*(zp3-zp2) elseif (b2c.lt.0.d0) then b2c=0.d0 b3c = (xp-xp1)*(xp3-xp1)+(yp-yp1)*(yp3-yp1)+(zp-zp1)*(zp3-zp1) b1c = (xp-xp3)*(xp1-xp3)+(yp-yp3)*(yp1-yp3)+(zp-zp3)*(zp1-zp3) endif b1c = max(b1c,0.d0) b2c = max(b2c,0.d0) b3c = max(b3c,0.d0) btc = b1c + b2c + b3c b1c = b1c / btc b2c = b2c / btc b3c = b3c / btc xn = b1c*xms(ip1)+b2c*xms(ip2)+b3c*xms(ip3) yn = b1c*yms(ip1)+b2c*yms(ip2)+b3c*yms(ip3) zn = b1c*zms(ip1)+b2c*zms(ip2)+b3c*zms(ip3) sn = sqrt (xn*xn + yn*yn + zn*zn) xn=xn/sn yn=yn/sn zn=zn/sn * verif compatibilite avec la normale au poin impactant if (xms(jp)*xn+yms(jp)*yn+zms(jp)*zn.gt.0.d0) goto 20 xjeu1 = (xp-xp1)*xn + (yp-yp1)*yn + (zp-zp1)*zn xjeu2 = (xp-xp2)*xn + (yp-yp2)*yn + (zp-zp2)*zn xjeu3 = (xp-xp3)*xn + (yp-yp3)*yn + (zp-zp3)*zn xjeu = xjeu1 * b1c + xjeu2 * b2c + xjeu3 * b3c ** write (6,*) ' normale ',xn,yn,zn,' jeu ',xjeu1,xjeu2,xjeu3 C*DBG write(ioimp,*) ' b1 b2 b3 ',b1,b2,b3 * creation du support du multiplicateur de lagrange si necessaire nbpts = xcoor(/1)/(idim+1) if (ipt1.num(1,iel).eq.ilgni) then nbpts = nbpts + 1 segadj mcoord segact ipt1*mod ipt1.num(1,iel)=nbpts endif segact mcoord*mod xcoor(nbpts*(idim+1)) = dist * Detection et traitement cas pathologiques if (b1.lt.dcrit) then if (b2.lt.dcrit) then C*DBG write(ioimp,*) 'path 1 (12)',ip3,jp nelpa1 = nelpa1 + 1 if (nelpa1.gt.nfopa1) then nfopa1 = nfopa1 + incnel segadj,mfopa1 endif lpoin(nelpa1,1)=ip3 lpoin(nelpa1,2)=jp lpoin(nelpa1,3)=ipt1.num(1,iel) xpoin(nelpa1,1)=xn xpoin(nelpa1,2)=yn xpoin(nelpa1,3)=zn else if (b3.lt.dcrit) then C*DBG write (ioimp,*) 'path 1 (13)',ip2,jp nelpa1 = nelpa1 + 1 if (nelpa1.gt.nfopa1) then nfopa1 = nfopa1 + incnel segadj,mfopa1 endif lpoin(nelpa1,1)=ip2 lpoin(nelpa1,2)=jp lpoin(nelpa1,3)=ipt1.num(1,iel) xpoin(nelpa1,1)=xn xpoin(nelpa1,2)=yn xpoin(nelpa1,3)=zn else C*DBG write (ioimp,*) 'path 2 (1)',ip2,ip3,jp nelpa2 = nelpa2 + 1 if (nelpa2.gt.nfopa2) then nfopa2 = nfopa2 + incnel segadj,mfopa2 endif lsegm(nelpa2,1)=min(ip2,ip3) lsegm(nelpa2,2)=max(ip2,ip3) lsegm(nelpa2,3)=jp lsegm(nelpa2,4)=ipt1.num(1,iel) xsegm(nelpa2,1)=xn xsegm(nelpa2,2)=yn xsegm(nelpa2,3)=zn endif goto 20 else if (b2.lt.dcrit) then if (b3.lt.dcrit) then C*DBG write(ioimp,*) 'path 1 (12)',ip3,jp nelpa1 = nelpa1 + 1 if (nelpa1.gt.nfopa1) then nfopa1 = nfopa1 + incnel segadj,mfopa1 endif lpoin(nelpa1,1)=ip1 lpoin(nelpa1,2)=jp lpoin(nelpa1,3)=ipt1.num(1,iel) xpoin(nelpa1,1)=xn xpoin(nelpa1,2)=yn xpoin(nelpa1,3)=zn else C*DBG write (ioimp,*) 'path 2 (2)',ip1,ip3,jp nelpa2 = nelpa2 + 1 if (nelpa2.gt.nfopa2) then nfopa2 = nfopa2 + incnel segadj,mfopa2 endif lsegm(nelpa2,1)=min(ip1,ip3) lsegm(nelpa2,2)=max(ip1,ip3) lsegm(nelpa2,3)=jp lsegm(nelpa2,4)=ipt1.num(1,iel) xsegm(nelpa2,1)=xn xsegm(nelpa2,2)=yn xsegm(nelpa2,3)=zn endif goto 20 else if (b3.lt.dcrit) then C*DBG write (ioimp,*) 'path 2 (3)',ip2,ip3,jp nelpa2 = nelpa2 + 1 if (nelpa2.gt.nfopa2) then nfopa2 = nfopa2 + incnel segadj,mfopa2 endif lsegm(nelpa2,1)=min(ip1,ip2) lsegm(nelpa2,2)=max(ip1,ip2) lsegm(nelpa2,3)=jp lsegm(nelpa2,4)=ipt1.num(1,iel) xsegm(nelpa2,1)=xn xsegm(nelpa2,2)=yn xsegm(nelpa2,3)=zn goto 20 endif * * Ajout d'une relation noeud-triangle * nelri3 = nelri3 + 1 nelch = nelch + 1 * * on ajuste les differents segments si necesaire * if (nelri3.gt.nelrig) then nelrig = nelrig + incnel segadj,xmatri nbelem = nbelem + incnel nbnn = 5 segadj,meleme nbnn = 1 segadj,ipt7 n = n + incnel segadj,mpoval endif * * Mise a jour du meleme num(1,nelri3) = ipt1.num(1,iel) num(2,nelri3) = ip1 num(3,nelri3) = ip2 num(4,nelri3) = ip3 num(5,nelri3) = jp * * Mise a jour de xmatri * lambda re( 1,1,nelri3) = 0d0 * ip1 re( 2,1,nelri3) = xn * b1 re( 3,1,nelri3) = yn * b1 re( 4,1,nelri3) = zn * b1 * ip2 re( 5,1,nelri3) = xn * b2 re( 6,1,nelri3) = yn * b2 re( 7,1,nelri3) = zn * b2 * ip3 re( 8,1,nelri3) = xn * b3 re( 9,1,nelri3) = yn * b3 re(10,1,nelri3) = zn * b3 * jp re(11,1,nelri3) = -xn re(12,1,nelri3) = -yn re(13,1,nelri3) = -zn * on transpose do 30 ic = 2, nligrp re(1,ic,nelri3) = re(ic,1,nelri3) 30 continue * le reste est nul * * remplissage du champoint de depi (jeu) * ipt7.num(1,nelch) = ipt1.num(1,iel) vpocha(nelch,1) = xjeu C*DBG write(ioimp,*) ' jeu ',xjeu,ip1,ip2,p3,jp,ipt1.num(1,iel) * 20 CONTINUE * * Ajustement au plus juste puis desactivation des segments lies * a la rigidite du type 3 if (nelri3.ne.nelrig) then nelrig = nelri3 segadj,xmatri nbelem = nelri3 nbnn = 5 segadj,meleme endif segdes,meleme,xmatri * write(ioimp,*) ' nb relation type 3 ',nelri3 C*DBG if (nelri3.eq.0) irigel(6,nrigel)=0 * * Relation du type 2 : noeud-segment *----------------------------------- if (nelpa2.eq.0) goto 200 * * Il y a au maximum nelpa2 relations de type 2 ! * * Ajustement du chpoint de depi (jeu) : mpoval et ipt7 * n = nelch + nelpa2 segadj,mpoval nbnn = 1 nbelem = nelch + nelpa2 nbsous = 0 nbref = 0 segadj,ipt7 * * creation de la raideur de conditions de type 2 * nrigel = nrigel+1 segadj,mrigid coerig(nrigel)=1.d0 * * creation du meleme associe a la relation * nbnn = 4 nbelem = nelpa2 nbsous = 0 nbref = 0 segini,meleme itypel=22 irigel(1,nrigel)=meleme * * creation du descriptif commun a toutes les raideurs * nligrp=10 nligrd=nligrp segini descr lisinc(1)='LX' lisdua(1)='FLX' noelep(1)=1 noeled(1)=1 do 110 i = 2, 10, 3 lisinc(i )=modepl(1) lisinc(i+1)=modepl(2) lisinc(i+2)=modepl(3) lisdua(i )=moforc(1) lisdua(i+1)=moforc(2) lisdua(i+2)=moforc(3) noelep(i )=(i+5)/3 noelep(i+1)=noelep(i) noelep(i+2)=noelep(i) noeled(i )=noelep(i) noeled(i+1)=noelep(i) noeled(i+2)=noelep(i) 110 continue segdes,descr irigel(3,nrigel)=descr * * creation du segment xmatri * nelrig=nelpa2 segini xmatri irigel(4,nrigel)=xmatri * * ce qu'on cree est unilateral * irigel(6,nrigel)=1 * * ce qu'on cree est symetrique * irigel(7,nrigel)=0 * * Nombre d'elements dans la rigidite du type 2 nelri2 = 0 * DO 120 iel = 1, nelpa2 * ip1 = lsegm(iel,1) if (ip1.eq.0) goto 120 C*DBG write (ioimp,*) 'point en relation ligne point ',ip1 lsegm(iel,1) = 0 ip2 = lsegm(iel,2) jp = lsegm(iel,3) il1 = lsegm(iel,4) xn = xsegm(iel,1) yn = xsegm(iel,2) zn = xsegm(iel,3) sn = sqrt(xn*xn+yn*yn+zn*zn) * DO 130 jel = iel+1,nelpa2 **pv goto 130 if (lsegm(jel,1).ne.ip1 .or. lsegm(jel,2).ne.ip2 .or. & lsegm(jel,3).ne.jp) goto 130 C*DBG write(ioimp,*) 'point en double ' * verif memes directions xn1 = xsegm(jel,1) yn1 = xsegm(jel,2) zn1 = xsegm(jel,3) if (xn*xn1+yn*yn1+zn*zn1 .lt. 0.9d0*sn) goto 130 lsegm(jel,1)=0 xn = xn + xn1 yn = yn + yn1 zn = zn + zn1 sn = sqrt(xn*xn+yn*yn+zn*zn) 130 CONTINUE if (.not.(sn.lt.1d-10) .and. .not.(sn.gt.1d-10)) & write(ioimp,*) 'Prob 6.1 - impo32' if (abs(sn).le.xpetit) write(ioimp,*) 'Prob 6.2 - impo32' C*DBG write(ioimp,*) ' xn yn zn sn ',xn,yn,zn,sn xn = xn / sn yn = yn / sn zn = zn / sn if (abs(xn).lt.1d-12) xn = 0.d0 if (abs(yn).lt.1d-12) yn = 0.d0 if (abs(zn).lt.1d-12) zn = 0.d0 * ipv = (ip1-1)*idimp1 xp1 = xcoor(ipv+1) yp1 = xcoor(ipv+2) zp1 = xcoor(ipv+3) ipv = (ip2-1)*idimp1 xp2 = xcoor(ipv+1) yp2 = xcoor(ipv+2) zp2 = xcoor(ipv+3) ipv = (jp -1)*idimp1 xp = xcoor(ipv+1) yp = xcoor(ipv+2) zp = xcoor(ipv+3) b1 = (xp-xp2)*(xp1-xp2)+(yp-yp2)*(yp1-yp2)+(zp-zp2)*(zp1-zp2) b2 = (xp-xp1)*(xp2-xp1)+(yp-yp1)*(yp2-yp1)+(zp-zp1)*(zp2-zp1) b1 = max(b1,0.d0) b2 = max(b2,0.d0) bt = b1+b2 if (.not.(bt.lt.1d-10) .and. .not.(bt.gt.1d-10)) & write(ioimp,*) 'Prob 7.1 - impo32' if (abs(bt).le.xpetit) write(ioimp,*) 'Prob 7.2 - impo32' b1 = b1 / bt b2 = b2 / bt * recalcul de la direction et du jeu en fonction des normales aux sommets xn = b1 *xms(ip1)+b2 *xms(ip2) yn = b1 *yms(ip1)+b2 *yms(ip2) zn = b1 *zms(ip1)+b2 *zms(ip2) sn = sqrt (xn*xn + yn*yn + zn*zn) xn=xn/sn yn=yn/sn zn=zn/sn xjeu1 = (xp-xp1)*xn + (yp-yp1)*yn + (zp-zp1)*zn xjeu2 = (xp-xp2)*xn + (yp-yp2)*yn + (zp-zp2)*zn xjeu = b1 * xjeu1 + b2 * xjeu2 C*DBG write(ioimp,*) 'type 2 ligne ',b1,b2 * nelri2 = nelri2 + 1 nelch = nelch + 1 * * on range dans le meleme * num(1,nelri2) = il1 num(2,nelri2) = ip1 num(3,nelri2) = ip2 num(4,nelri2) = jp * * on remplit le xmatri * * lambda re( 1,1,nelri2) = 0d0 * ip1 re( 2,1,nelri2) = xn*b1 re( 3,1,nelri2) = yn*b1 re( 4,1,nelri2) = zn*b1 * ip2 re( 5,1,nelri2) = xn*b2 re( 6,1,nelri2) = yn*b2 re( 7,1,nelri2) = zn*b2 * jp re( 8,1,nelri2) = -xn re( 9,1,nelri2) = -yn re(10,1,nelri2) = -zn * on transpose do 140 ic=2,10 re(1,ic,nelri2) = re(ic,1,nelri2) 140 continue * le reste est nul * * remplissage du champoint de depi (jeu) * ipt7.num(1,nelch) = il1 vpocha(nelch,1) = xjeu C*DBG write(ioimp,*) ' jeu ',xjeu,ip1,ip2,jp,il1,il2 * 120 CONTINUE * * Ajustement au plus juste puis desactivation des segments lies * a la rigidite du type 2 if (nelri2.ne.nelrig) then nelrig = nelri2 segadj xmatri nbelem = nelri2 nbnn = 4 segadj,meleme endif segdes,xmatri,meleme * write (ioimp,*) ' nb relation type 2 ',nelri2 C*DBG if (nelri2.eq.0) irigel(6,nrigel)=0 * * Relation du type 1 : noeud-noeud *---------------------------------- 200 CONTINUE if (nelpa1.eq.0) goto 300 * * Il y a au maximum nelpa1 relations de type 1 ! * * Ajustement du chpoint de depi (jeu) : mpoval et ipt7 * n = nelch + nelpa1 segadj,mpoval nbnn = 1 nbelem = nelch + nelpa1 nbsous = 0 nbref = 0 segadj,ipt7 * * creation de la raideur de conditions * nrigel = nrigel + 1 segadj mrigid coerig(nrigel) = 1.d0 * * creation du meleme associe a la relation * nbnn = 3 nbelem = nelpa1 nbsous = 0 nbref = 0 segini meleme itypel = 22 irigel(1,nrigel) = meleme * * creation du descriptif commun a toutes les raideurs * nligrp=7 nligrd=nligrp segini descr lisinc(1)='LX ' lisdua(1)='FLX ' noelep(1)=1 noeled(1)=1 do 210 i=2,7,3 lisinc(i )=modepl(1) lisinc(i+1)=modepl(2) lisinc(i+2)=modepl(3) lisdua(i )=moforc(1) lisdua(i+1)=moforc(2) lisdua(i+2)=moforc(3) noelep(i )=(i+5)/3 noelep(i+1)=noelep(i) noelep(i+2)=noelep(i) noeled(i )=noelep(i) noeled(i+1)=noelep(i) noeled(i+2)=noelep(i) 210 continue segdes,descr irigel(3,nrigel) = descr * * creation du segment xmatri * nelrig = nelpa1 segini,xmatri irigel(4,nrigel)=xmatri * * ce qu'on cree est unilateral * irigel(6,nrigel)=1 * * ce qu'on cree est symetrique * irigel(7,nrigel)=0 * * Nombre d'elements dans la rigidite du type 1 nelri1 = 0 * DO 220 iel = 1, nelpa1 * ip1 = lpoin(iel,1) if (ip1.eq.0) goto 220 lpoin(iel,1) = 0 jp = lpoin(iel,2) il1 = lpoin(iel,3) * xn = xpoin(iel,1) yn = xpoin(iel,2) zn = xpoin(iel,3) sn = sqrt(xn*xn+yn*yn+zn*zn) C*DBG write (ioimp,*) 'point en relation point point ', C*DBG & ip1,jp,il1,xn,yn,zn DO 230 jel = iel+1, nelpa1 if (lpoin(jel,1).ne.ip1.or.lpoin(jel,2).ne.jp) goto 231 * verif memes directions xn1 = xpoin(jel,1) yn1 = xpoin(jel,2) zn1 = xpoin(jel,3) **pv if (xn*xn1+yn*yn1+zn*zn1 .lt. 0.9d0*sn) goto 230 C*DBG write (ioimp,*) 'point en double cas 1',ip1,jp lpoin(jel,1) = 0 xn = xn + xn1 yn = yn + yn1 zn = zn + zn1 sn = sqrt(xn*xn+yn*yn+zn*zn) 231 continue **pv goto 230 if (lpoin(jel,1).ne.jp.or.lpoin(jel,2).ne.ip1) goto 230 * verif memes directions xn1 = xpoin(jel,1) yn1 = xpoin(jel,2) zn1 = xpoin(jel,3) if (xn*xn1+yn*yn1+zn*zn1 .gt. -0.9d0*sn) goto 230 C*DBG write (ioimp,*) 'point en double cas 2',jp,ip1 lpoin(jel,1) = 0 xn = xn - xn1 yn = yn - yn1 zn = zn - zn1 sn = sqrt(xn*xn+yn*yn+zn*zn) 230 CONTINUE C*DBG write (ioimp,*) ' xn yn zn sn ',xn,yn,zn,sn if (.not.(sn.lt.1d-10).and. .not.(sn.gt.1d-10)) & write(ioimp,*) 'Prob 8.1 - impo32' if (abs(sn).le.xpetit) write(ioimp,*) 'Prob 8.2 - impo32' xn = xn / sn yn = yn / sn zn = zn / sn if (abs(xn).lt.1d-12) xn = 0.d0 if (abs(yn).lt.1d-12) yn = 0.d0 if (abs(zn).lt.1d-12) zn = 0.d0 * ipv = (ip1-1)*idimp1 xp1 = xcoor(ipv+1) yp1 = xcoor(ipv+2) zp1 = xcoor(ipv+3) ipv = (jp -1)*idimp1 xp = xcoor(ipv+1) yp = xcoor(ipv+2) zp = xcoor(ipv+3) * recalcul de la direction et du jeu en fonction des normales aux sommets xn = xms(ip1) yn = yms(ip1) zn = zms(ip1) sn = sqrt (xn*xn + yn*yn + zn*zn) xn=xn/sn yn=yn/sn zn=zn/sn xjeu = (xp-xp1)*xn + (yp-yp1)*yn + (zp-zp1)*zn * nelri1 = nelri1 + 1 nelch = nelch + 1 * * on range dans le meleme * num(1,nelri1) = il1 num(2,nelri1) = ip1 num(3,nelri1) = jp * * on remplit le xmatri * * lambda re(1,1,nelri1) = 0d0 * ip1 re(2,1,nelri1) = xn re(3,1,nelri1) = yn re(4,1,nelri1) = zn * jp re(5,1,nelri1) = -xn re(6,1,nelri1) = -yn re(7,1,nelri1) = -zn * on transpose do 240 ic=2,7 re(1,ic,nelri1) = re(ic,1,nelri1) 240 continue * le reste est nul * * remplissage du champoint de depi (jeu) ipt7.num(1,nelch) = il1 vpocha(nelch,1) = xjeu C*DBG write(ioimp,*) ' jeu ',xjeu,ip1,jp,il1,il2 * 220 CONTINUE * * Ajustement au plus juste puis desactivation des segments lies * a la rigidite du type 1 if (nelri1.ne.nelrig) then nelrig = nelri1 segadj xmatri nbelem = nelri1 nbnn = 3 segadj,meleme endif segdes,xmatri,meleme * C*DBG if (nelri1.eq.0) irigel(6,nrigel)=0 write(ioimp,*) ' nb relation type 1 ',nelri1 * * Fin du traitement de la formulation standard/forte * ---------------------------------------------------- 300 CONTINUE * Destruction des segments locaux devenus inutiles segsup,mfopa1,mfopa2,mfopa3 GOTO 1000 *======================================================================= *= Formulation "faible" 3D : relation triangle-triangle *======================================================================= * Element de contact 3D a 7 noeuds (1+6) 500 CONTINUE * * Petit segment de travail segini,mfaible * * Relation du type 0 : triangle-triangle (faible) *--------------------------------------- * * Creation du meleme associe a la relation * nbnn = 7 nbelem = incnel nbsous = 0 nbref = 0 segini,meleme itypel = 22 irigel(1,nrigel) = meleme * * Creation du descriptif commun a toutes les raideurs * nligrp = 19 nligrd = nligrp segini,descr lisinc(1) = 'LX ' lisdua(1) = 'FLX ' noelep(1) = 1 noeled(1) = 1 do 510 i = 2, nligrp, 3 lisinc(i ) = modepl(1) lisinc(i+1) = modepl(2) lisinc(i+2) = modepl(3) lisdua(i ) = moforc(1) lisdua(i+1) = moforc(2) lisdua(i+2) = moforc(3) noelep(i ) = (i+5)/3 noelep(i+1) = noelep(i) noelep(i+2) = noelep(i) noeled(i ) = noelep(i) noeled(i+1) = noelep(i) noeled(i+2) = noelep(i) 510 continue segdes,descr irigel(3,nrigel) = descr * * creation du segment xmatri * nelrig = incnel segini,xmatri irigel(4,nrigel) = xmatri * * ce qu'on cree est unilateral * irigel(6,nrigel) = 1 * * ce qu'on cree est symetrique * irigel(7,nrigel) = 0 * * Nombre d'elements crees dans meleme=irigel(1,nrigel), ipt7 et mpoval nelri0 = 0 * * Boucle sur les elements du maillage de contact/frottement "faible" * DO 520 iel = 1, nbel1 * ip1 = ipt1.num(2,iel) ip2 = ipt1.num(3,iel) ip3 = ipt1.num(4,iel) jp1 = ipt1.num(5,iel) jp2 = ipt1.num(6,iel) jp3 = ipt1.num(7,iel) C*DBG write(ioimp,*) iel,ip1,ip2,ip3,jp1,jp2,jp3,ipt1.num(1,iel) * * Verification (provisoire) qu'il n'y a pas de noeud double if (ip1.eq.jp1) goto 520 if (ip1.eq.jp2) goto 520 if (ip1.eq.jp3) goto 520 if (ip2.eq.jp1) goto 520 if (ip2.eq.jp2) goto 520 if (ip2.eq.jp3) goto 520 if (ip3.eq.jp1) goto 520 if (ip3.eq.jp2) goto 520 if (ip3.eq.jp3) goto 520 * * Definition du triangle T1(ip1,ip2,ip3) * - Coordonnees des noeuds de T1 ipv = (ip1-1)*idimp1 xp1 = xcoor(ipv+1) yp1 = xcoor(ipv+2) zp1 = xcoor(ipv+3) ipv = (ip2-1)*idimp1 xp2 = xcoor(ipv+1) yp2 = xcoor(ipv+2) zp2 = xcoor(ipv+3) ipv = (ip3-1)*idimp1 xp3 = xcoor(ipv+1) yp3 = xcoor(ipv+2) zp3 = xcoor(ipv+3) * - Barycentre de T1 xgT1 = (xp1 + xp2 + xp3) * X1s3 ygT1 = (yp1 + yp2 + yp3) * X1s3 zgT1 = (zp1 + zp2 + zp3) * X1s3 * - Normale au triangle T1 xnT1 = (yp1-yp2)*(zp2-zp3) - (zp1-zp2)*(yp2-yp3) ynT1 = (zp1-zp2)*(xp2-xp3) - (xp1-xp2)*(zp2-zp3) znT1 = (xp1-xp2)*(yp2-yp3) - (yp1-yp2)*(xp2-xp3) dnT1 = SQRT(xnT1*xnT1+ynT1*ynT1+znT1*znT1) IF (.not.(dnT1.lt.1d-10) .and. .not.(dnT1.gt.1d-10)) & write(ioimp,*) 'FAIBle - 1.1 - impo32' IF (abs(dnT1).le.xpetit) write(ioimp,*) 'FAIBle - 1.2 - impo32' xnT1 = xnT1 / dnT1 ynT1 = ynT1 / dnT1 znT1 = znT1 / dnT1 if (abs(xnT1).lt.1d-12) xnT1 = 0.d0 if (abs(ynT1).lt.1d-12) ynT1 = 0.d0 if (abs(znT1).lt.1d-12) znT1 = 0.d0 C*DBG write(ioimp,*) 'Normale au plan T1',xnT1,ynT1,znT1,dnT1 * - Rayon de la sphere "englobante" de l'element (centre=barycentre) d1 = (xgT1-xp1)**2 + (ygT1-yp1)**2 + (zgT1-zp1)**2 d2 = (xgT1-xp2)**2 + (ygT1-yp2)**2 + (zgT1-zp2)**2 d3 = (xgT1-xp3)**2 + (ygT1-yp3)**2 + (zgT1-zp3)**2 RayT1 = SQRT(max(d1,d2,d3)) * * Definition du triangle T2(jp1,jp2,jp3) * - Coordonnees des noeuds de T2 ipv = (jp1-1)*idimp1 xq1 = xcoor(ipv+1) yq1 = xcoor(ipv+2) zq1 = xcoor(ipv+3) ipv = (jp2-1)*idimp1 xq2 = xcoor(ipv+1) yq2 = xcoor(ipv+2) zq2 = xcoor(ipv+3) ipv = (jp3-1)*idimp1 xq3 = xcoor(ipv+1) yq3 = xcoor(ipv+2) zq3 = xcoor(ipv+3) * - Barycentre de T2 xgT2 = (xq1 + xq2 + xq3) * X1s3 ygT2 = (yq1 + yq2 + yq3) * X1s3 zgT2 = (zq1 + zq2 + zq3) * X1s3 * - Normale au triangle T2 xnT2 = (yq1-yq2)*(zq2-zq3) - (zq1-zq2)*(yq2-yq3) ynT2 = (zq1-zq2)*(xq2-xq3) - (xq1-xq2)*(zq2-zq3) znT2 = (xq1-xq2)*(yq2-yq3) - (yq1-yq2)*(xq2-xq3) dnT2 = SQRT(xnT2*xnT2+ynT2*ynT2+znT2*znT2) IF (.not.(dnT2.lt.1d-10) .and. .not.(dnT2.gt.1d-10)) & write(ioimp,*) 'FAIBle - 2.1 - impo32' IF (abs(dnT2).le.xpetit) write(ioimp,*) 'FAIBle - 2.2 - impo32' xnT2 = xnT2 / dnT2 ynT2 = ynT2 / dnT2 znT2 = znT2 / dnT2 if (abs(xnT2).lt.1d-12) xnT2 = 0.d0 if (abs(ynT2).lt.1d-12) ynT2 = 0.d0 if (abs(znT2).lt.1d-12) znT2 = 0.d0 C*DBG write(ioimp,*) 'Normale au plan T2',xnT2,ynT2,znT2,dnT2 * - Rayon de la sphere "englobante" de l'element (centre=barycentre) d1 = (xgT2-xq1)**2 + (ygT2-yq1)**2 + (zgT2-zq1)**2 d2 = (xgT2-xq2)**2 + (ygT2-yq2)**2 + (zgT2-zq2)**2 d3 = (xgT2-xq3)**2 + (ygT2-yq3)**2 + (zgT2-zq3)**2 RayT2 = SQRT(max(d1,d2,d3)) * Orientation respective correcte des 2 normales * Les triangles T1 et T2 etant decrits dans le meme sens ("prisme"), * le produit scalaire de leur normale doit etre positif ou nul ! scal = xnT1 * xnT2 + ynT1 * ynT2 + znT1 * znT2 C*DBG write(ioimp,*) 'Prod.scal des normales',scal if (scal.lt.0.) goto 520 * * Distance entre les centres de gravite * Les spheres "englobantes" doivent s'intersecter pour contact potentiel dist = SQRT((xgT2-xgT1)**2 + (ygT2-ygT1)**2 + (zgT2-zgT1)**2) IF (dist.GT.RayT1+RayT2) goto 520 * * Critere d'acceptation de l'element de contact * * Definition du plan de contact (surface intermediaire) * - Normale (unitaire) du plan de contact xnC = xnT1 + xnT2 ynC = ynT1 + ynT2 znC = znT1 + znT2 dnC = SQRT(xnC*xnC + ynC*ynC + znC*znC) IF (.not.(dnC.lt.1d-10) .and. .not.(dnC.GT.1D-10)) & write(ioimp,*) 'FAIBle - 3.1 - impo32' IF (abs(dnC).le.xpetit) write(ioimp,*) 'FAIBle - 3.2 - impo32' xnC = xnC / dnC ynC = ynC / dnC znC = znC / dnC if (abs(xnC).lt.1d-12) xnC = 0.d0 if (abs(ynC).lt.1d-12) ynC = 0.d0 if (abs(znC).lt.1d-12) znC = 0.d0 C*DBG write(ioimp,*) 'Normale au plan Contact',xnC,ynC,znC,dnC * - Point definissant le plan de contact * = milieu des barycentres des triangles T1 et T2 xgC = (xgT1 + xgT2) * 0.5 ygC = (ygT1 + ygT2) * 0.5 zgC = (zgT1 + zgT2) * 0.5 * - "Distance" de reference au plan de contact dgC = xgC*xnC + ygC*ynC + zgC*znC * * Distance "signee" des triangles T1 et T2 au plan de contact dp1C = (xp1*xnC + yp1*ynC + zp1*znC) - dgC dp2C = (xp2*xnC + yp2*ynC + zp2*znC) - dgC dp3C = (xp3*xnC + yp3*ynC + zp3*znC) - dgC dq1C = (xq1*xnC + yq1*ynC + zq1*znC) - dgC dq2C = (xq2*xnC + yq2*ynC + zq2*znC) - dgC dq3C = (xq3*xnC + yq3*ynC + zq3*znC) - dgC * Projection des triangles T1 et T2 sur le plan de contact * T1C = triangle T1 projete sur plan de Contact cPT1C(1,1) = xp1 - dp1C * xnC cPT1C(1,2) = yp1 - dp1C * ynC cPT1C(1,3) = zp1 - dp1C * znC cPT1C(2,1) = xp2 - dp2C * xnC cPT1C(2,2) = yp2 - dp2C * ynC cPT1C(2,3) = zp2 - dp2C * znC cPT1C(3,1) = xp3 - dp3C * xnC cPT1C(3,2) = yp3 - dp3C * ynC cPT1C(3,3) = zp3 - dp3C * znC * T2C = triangle T2 projete sur plan de Contact cPT2C(1,1) = xq1 - dq1C * xnC cPT2C(1,2) = yq1 - dq1C * ynC cPT2C(1,3) = zq1 - dq1C * znC cPT2C(2,1) = xq2 - dq2C * xnC cPT2C(2,2) = yq2 - dq2C * ynC cPT2C(2,3) = zq2 - dq2C * znC cPT2C(3,1) = xq3 - dq3C * xnC cPT2C(3,2) = yq3 - dq3C * ynC cPT2C(3,3) = zq3 - dq3C * znC * * On determine quelle est la composante maximale de la normale pour * projeter les triangles T1C et T2C sur un plan parallele aux axes * ("le plus orthogonal" a cette normale), ce qui maximise leur aire. r_x = abs(xnC) r_y = abs(ynC) r_z = abs(znC) if (r_x.gt.r_y) then inC = 1 if (r_x.lt.r_z) inC = 3 else inC = 2 if (r_y.lt.r_z) inC = 3 endif * Projection des triangles T1C et T2C sur ce plan // aux axes (A) * T1A = triangle T1C projete sur ce plan (A) * T2A = triangle T2C projete sur ce plan (A) * On prend soin d'orienter les triangles T1A et T2A dans le sens * direct (en tenant compte du signe de la composante normale !) GOTO (531,532,533), inC call erreur(5) return * Normale NC selon "x" : A =(y,z) 531 if (xnC.lt.0.) inC = -inC inX = 2 inY = 3 goto 534 * Normale NC selon "y" : A =(x,z) 532 if (ynC.lt.0.) inC = -inC inX = 1 inY = 3 goto 534 * Normale NC selon "z" : A =(x,y) 533 if (znC.lt.0.) inC = -inC inX = 1 inY = 2 goto 534 534 continue cPT1A(1,1) = cPT1C(1,inX) cPT1A(1,2) = cPT1C(1,inY) cPT2A(1,1) = cPT2C(1,inX) cPT2A(1,2) = cPT2C(1,inY) if (inC.gt.0) then cPT1A(2,1) = cPT1C(2,inX) cPT1A(2,2) = cPT1C(2,inY) cPT2A(2,1) = cPT2C(2,inX) cPT2A(2,2) = cPT2C(2,inY) cPT1A(3,1) = cPT1C(3,inX) cPT1A(3,2) = cPT1C(3,inY) cPT2A(3,1) = cPT2C(3,inX) cPT2A(3,2) = cPT2C(3,inY) else cPT1A(2,1) = cPT1C(3,inX) cPT1A(2,2) = cPT1C(3,inY) cPT2A(2,1) = cPT2C(3,inX) cPT2A(2,2) = cPT2C(3,inY) cPT1A(3,1) = cPT1C(2,inX) cPT1A(3,2) = cPT1C(2,inY) cPT2A(3,1) = cPT2C(2,inX) cPT2A(3,2) = cPT2C(2,inY) endif cPT1A(4,1) = cPT1A(1,1) cPT1A(4,2) = cPT1A(1,2) * Calcul des cotes des triangles T1A et T2A (12,23,31) * Utile pour connaitre ensuite la normale aux cotes des triangles vT1A(1,1) = cPT1A(2,1) - cPT1A(1,1) vT1A(1,2) = cPT1A(2,2) - cPT1A(1,2) vT1A(2,1) = cPT1A(3,1) - cPT1A(2,1) vT1A(2,2) = cPT1A(3,2) - cPT1A(2,2) vT1A(3,1) = cPT1A(1,1) - cPT1A(3,1) vT1A(3,2) = cPT1A(1,2) - cPT1A(3,2) C*nu vT2A(1,1) = cPT2A(2,1) - cPT2A(1,1) vT2A(1,2) = cPT2A(2,2) - cPT2A(1,2) C*nu vT2A(2,1) = cPT2A(3,1) - cPT2A(2,1) vT2A(2,2) = cPT2A(3,2) - cPT2A(2,2) C*nu vT2A(3,1) = cPT2A(1,1) - cPT2A(3,1) vT2A(3,2) = cPT2A(1,2) - cPT2A(3,2) * Calcul de la surface de chacun des triangles T1A et T2A * (en fait, on calcule le double, mais par la suite on ne considere que * des rapports de surfaces de triangle) SuT1A = (cPT1A(2,1)+cPT1A(1,1)) * vT1A(1,2) & + (cPT1A(3,1)+cPT1A(2,1)) * vT1A(2,2) & + (cPT1A(1,1)+cPT1A(3,1)) * vT1A(3,2) SuT2A = (cPT2A(2,1)+cPT2A(1,1)) * vT2A(1,2) & + (cPT2A(3,1)+cPT2A(2,1)) * vT2A(2,2) & + (cPT2A(1,1)+cPT2A(3,1)) * vT2A(3,2) C*DBG write(ioimp,*) 'Surfaces T1A - T2A :', 0.5*SuT1A,0.5*SuT2A C*DBG if (SuT1A.gt.100.*SuT2A .or. SuT2A.gt.100.*SuT1A) C*DBG & write(ioimp,*) 'Rapport des surfaces tres important !' * On initialise l'intersection des triangles avec T2A nPIn0 = 3 do 540 i = 1, nPIn0 cPIn0(i,1) = cPT2A(i,1) cPIn0(i,2) = cPT2A(i,2) 540 continue * Determination de l'intersection des triangles T1A et T2A * en regardant progressivement l'intersection avec les cotes de T1A * Note : On utilise le fait que les polygones traites sont convexes ! DO 550 i = 1, 3 vN1x = -vT1A(i,2) vN1y = vT1A(i,1) scal = vN1x * cPT1A(i,1) + vN1y * cPT1A(i,2) * ipos = 0 ineg = 0 jpos = 0 jneg = 0 do 551 j = 1, nPIn0 test(j) = (vN1x * cPIn0(j,1) + vN1y * cPIn0(j,2)) - scal if (test(j).gt.1.d-10) then * if (test(j).gt.0.) then ipos = ipos + 1 if (jneg.ne.0 .and. jpos.eq.0) jpos = j else if (test(j).lt.-1.d-10) then * else if (test(j).lt.0.) then ineg = ineg + 1 if (jneg.eq.0) jneg = j else test(j) = 0. endif 551 continue * * 1) Cas ou ipos = 0 : * 1.1) il n'y pas d'intersection (ineg=nPIn0) * 1.2) l'intersection se limite a 1 point (ineg=nPIn0-1) * 1.3) l'intersection a 1 segment sur une frontiere (ineg=nPIn0-2) * La surface correspondante est nulle et il n'y a donc pas de contact ! if (ipos.eq.0) then C*DBG write(ioimp,*) 'Intersection a ',nPIn0-ineg,' points' goto 520 endif * 2) Cas ou ipos > 0 : Il y a une intersection a calculer * 2.1) ineg = 0 : Tous les points sont "du bon cote" et l'intersection * correspond au polygone de depart ! if (ineg.eq.0) then C*DBG write(ioimp,*) 'Intersection du bon cote - rien a faire' goto 550 endif * 2.2) ineg > 0 : Il y a deux intersections a determiner car il y a * au moins un point du "mauvais cote" jpos = max(1,jpos) jpr = jpos - 1 if (jpos.eq.1) jpr = nPIn0 nPIn = 1 if (test(jpr).lt.0.) then r_z = test(jpos) / (test(jpos) - test(jpr)) cPIn(nPIn,1) = cPIn0(jpos,1) & + r_z*(cPIn0(jpr,1) - cPIn0(jpos,1)) cPIn(nPIn,2) = cPIn0(jpos,2) & + r_z*(cPIn0(jpr,2) - cPIn0(jpos,2)) else cPIn(nPIn,1) = cPIn0(jpr,1) cPIn(nPIn,2) = cPIn0(jpr,2) endif do 552 j = 1, ipos nPIn = nPIn + 1 cPIn(nPIn,1) = cPIn0(jpos,1) cPIn(nPIn,2) = cPIn0(jpos,2) jpr = jpos jpos = jpos + 1 if (jpos.gt.nPIn0) jpos = 1 552 continue nPIn = nPIn + 1 if (test(jpos).lt.0.) then r_z = test(jpr) / (test(jpr) - test(jpos)) cPIn(nPIn,1) = cPIn0(jpr,1) & + r_z*(cPIn0(jpos,1) - cPIn0(jpr,1)) cPIn(nPIn,2) = cPIn0(jpr,2) & + r_z*(cPIn0(jpos,2) - cPIn0(jpr,2)) else cPIn(nPIn,1) = cPIn0(jpos,1) cPIn(nPIn,2) = cPIn0(jpos,2) endif * Mise a jour cPIn0 pour la suite du traitement des intersections... nPIn0 = nPIn do 554 j = 1, nPIn0 cPIn0(j,1) = cPIn(j,1) cPIn0(j,2) = cPIn(j,2) 554 continue C*DBG-F write(ioimp,*) 'Intersection ',i,' a ',nPIn0,' points' C*DBG-F do j = 1, nPIn0 C*DBG-F write(ioimp,*) ' ',j,nbpts0+6+j,cPIn0(j,1),cPIn0(j,2) C*DBG-F enddo * 550 CONTINUE * * Calcul de la surface de l'intersection et de son centre de gravite r_z = cPIn0(nPIn0,1)*cPIn0(1,2) - cPIn0(1,1)*cPIn0(nPIn0,2) SuPIn = r_z xGIn = (cPIn0(nPIn0,1)+cPIn0(1,1)) * r_z yGIn = (cPIn0(nPIn0,2)+cPIn0(1,2)) * r_z do 560 i = 1, nPIn0-1 j = i + 1 r_z = cPIn0(i,1)*cPIn0(j,2) - cPIn0(j,1)*cPIn0(i,2) SuPIn = SuPIn + r_z xGIn = xGIn + (cPIn0(i,1)+cPIn0(j,1)) * r_z yGIn = yGIn + (cPIn0(i,2)+cPIn0(j,2)) * r_z 560 continue if (SuPIn .lt. 1.E-4*max(SuT1A,SuT2A) ) then C*DBG write(ioimp,*) 'Intersection a une surface negligeable !' goto 520 endif r_z = X1s3 / SupIn SuPIn = 0.5 * SupIn xGIn = xGIn * r_z yGIn = yGIn * r_z C*DBG write(ioimp,*) 'Surface Intersection :',SuPIn * * Calcul des coordonnees barycentriques du centre de gravite * pour le triangle T1A xpip1 = xGIn - cPT1A(1,1) ypip1 = yGIn - cPT1A(1,2) xpip2 = xGIn - cPT1A(2,1) ypip2 = yGIn - cPT1A(2,2) xpip3 = xGIn - cPT1A(3,1) ypip3 = yGIn - cPT1A(3,2) b1T1 = xpip2 * ypip3 - ypip2 * xpip3 b2T1 = xpip3 * ypip1 - ypip3 * xpip1 b3T1 = xpip1 * ypip2 - ypip1 * xpip2 bt = b1T1 + b2T1 + b3T1 if (abs(bt-SuT1A) .gt. 1.E-3) & write(ioimp,*) 'Prob bt-SuT1A',bt,SuT1A b1T1 = b1T1 / SuT1A b2T1 = b2T1 / SuT1A b3T1 = b3T1 / SuT1A * * pour le triangle T2A xpip1 = xGIn - cPT2A(1,1) ypip1 = yGIn - cPT2A(1,2) xpip2 = xGIn - cPT2A(2,1) ypip2 = yGIn - cPT2A(2,2) xpip3 = xGIn - cPT2A(3,1) ypip3 = yGIn - cPT2A(3,2) b1T2 = xpip2 * ypip3 - ypip2 * xpip3 b2T2 = xpip3 * ypip1 - ypip3 * xpip1 b3T2 = xpip1 * ypip2 - ypip1 * xpip2 bt = b1T2 + b2T2 + b3T2 if (abs(bt-SuT2A) .gt. 1.E-3) & write(ioimp,*) 'Prob bt-SuT2A',bt,SuT2A b1T2 = b1T2 / SuT2A b2T2 = b2T2 / SuT2A b3T2 = b3T2 / SuT2A * * Calcul du jeu xjeu = (b1T2*dq1C + b2T2*dq2C + b3T2*dq3C) & - (b1T1*dp1C + b2T1*dp2C + b3T1*dp3C) C*DBG write(ioimp,*) 'Jeu =',xjeu * * Ajout d'une relation triangle-triangle * nelri0 = nelri0 + 1 nelch = nelch + 1 * * on ajuste les differents segments si necesaire * if (nelri0.gt.nelrig) then nelrig = nelrig + incnel segadj,xmatri nbelem = nbelem + incnel nbnn = 7 segadj,meleme nbnn = 1 segadj,ipt7 n = n + incnel segadj,mpoval endif * * Mise a jour du meleme * if (ipt1.num(1,iel).eq.ILGNI) then segact ipt1*mod nbpts=xcoor(/1)/(idim+1)+1 segadj mcoord ipt1.num(1,iel)=nbpts endif num(1,nelri0) = ipt1.num(1,iel) num(2,nelri0) = ip1 num(3,nelri0) = ip2 num(4,nelri0) = ip3 num(5,nelri0) = jp1 num(6,nelri0) = jp2 num(7,nelri0) = jp3 *$ * Mise a jour de xmatri * * lambda re( 1,1,nelri0) = 0.d0 * ip1 re( 2,1,nelri0) = xnC * b1T1 re( 3,1,nelri0) = ynC * b1T1 re( 4,1,nelri0) = znC * b1T1 * ip2 re( 5,1,nelri0) = xnC * b2T1 re( 6,1,nelri0) = ynC * b2T1 re( 7,1,nelri0) = znC * b2T1 * ip3 re( 8,1,nelri0) = xnC * b3T1 re( 9,1,nelri0) = ynC * b3T1 re(10,1,nelri0) = znC * b3T1 * jp1 re(11,1,nelri0) = -xnC * b1T2 re(12,1,nelri0) = -ynC * b1T2 re(13,1,nelri0) = -znC * b1T2 * jp2 re(14,1,nelri0) = -xnC * b2T2 re(15,1,nelri0) = -ynC * b2T2 re(16,1,nelri0) = -znC * b2T2 * jp3 re(17,1,nelri0) = -xnC * b3T2 re(18,1,nelri0) = -ynC * b3T2 re(19,1,nelri0) = -znC * b3T2 * on transpose do 580 ic = 2, nligrp re(1,ic,nelri0) = re(ic,1,nelri0) 580 continue * le reste est nul * * remplissage du champoint de depi (jeu) * ipt7.num(1,nelch) = ipt1.num(1,iel) vpocha(nelch,1) = xjeu C*DBG-F call prfaible(ipt1,iel,nPIn0,mfaible) * 520 CONTINUE * * Ajustement au plus juste puis desactivation des segments lies * a la rigidite du type 0 if (nelri0.ne.nelrig) then nelrig = nelri0 segadj,xmatri nbelem = nelri0 nbnn = 7 segadj,meleme endif segdes,meleme,xmatri * write(ioimp,*) ' nb relation type 0 ',nelri0 * if (nelri0.eq.0) irigel(6,nrigel)=0 * * Destruction des segments locaux devenus inutiles segsup,mfaible * C* GOTO 1000 *======================================================================= *= Fin du sous-programme : *======================================================================= 1000 CONTINUE * * Ajustement au plus juste du chpoint de depi (jeu) : mpoval et ipt7 * puis desactivation du chpoint if (n.ne.nelch) then n = nelch segadj,mpoval nbnn = 1 nbelem = nelch nbsous = 0 nbref = 0 segadj,ipt7 endif segdes,mpoval,ipt7 segdes,msoupo,mchpoi * Desctivation de la matrice de raideur de contact segdes,mrigid * * Ecriture des resultats de l'operateur : call ecrobj('RIGIDITE',mrigid) call ecrobj('CHPOINT ',mchpoi) * 900 CONTINUE segdes,ipt1 return end