Télécharger la plaquette Cast3M


fr en

Télécharger impos2.eso

Retour à la liste

Numérotation des lignes : 

impos2

  1. C IMPOS2    SOURCE    MB234859  24/12/16    21:15:03     12099          
  2.  
  3. *  impo bloc en 2D
  4.  
  5.       SUBROUTINE IMPOS2(ipt1,ipt6,ipt8,itcont,mchel1,mrigid,mchpoi)
  6.  
  7.       IMPLICIT INTEGER(I-N)
  8.       IMPLICIT REAL*8 (A-H,O-Z)
  9.  
  10. -INC PPARAM
  11. -INC CCOPTIO
  12.  
  13. -INC SMCOORD
  14. -INC SMELEME
  15. -INC SMRIGID
  16. -INC SMCHPOI
  17. -INC SMCHAML
  18. -INC CCREEL
  19.  
  20. * directions moyennes aux sommets
  21.       segment mfopa2
  22.         real*8 xms(nbpts),yms(nbpts)
  23.       endsegment
  24.  
  25.       character*4 modepl(4),moforc(4)
  26.  
  27.       data modepl /'UX  ','UY  ','UR  ','UZ  '/
  28.       data moforc /'FX  ','FY  ','FR  ','FZ  '/
  29.  
  30. * definition de tableaux utilises pour mortar
  31.       real*8 xiG(2)
  32.       real*8 zetaG(2)
  33.       real*8 wg(2)
  34. *
  35.       idimp1 = IDIM + 1
  36. *
  37. *  Activation du maillage et petites verifications
  38. *
  39.       segact ipt1,ipt6,ipt8
  40.       nbno1 = ipt1.num(/1)
  41.       nbel1 = ipt1.num(/2)
  42.       nbno6 = ipt6.num(/1)
  43.       nbel6 = ipt6.num(/2)
  44.       nbno8 = ipt8.num(/1)
  45.       nbel8 = ipt8.num(/2)
  46.       if (ipt1.lisous(/1).ne.0) call erreur(25)
  47. **    write(6,*) 'nbno1 nbno6 nbno8 ',nbno1,nbno6,nbno8
  48.       if (ipt6.itypel.ne.2) call erreur(16)
  49.       if (nbno6.ne.2) call erreur(16)
  50.       if (ierr.ne.0) goto 900
  51. *
  52. *  Indice des ddls concernes en fonction du mode 2D
  53. *
  54.       imo = 1
  55.       if (ifour.eq.0) imo = 3
  56. *
  57. *   Increment sur le nombre d'elements de contact retenus et utilises
  58. *   dans le chpoint (mpoval et igeoc) et rigidite pour agrandir les
  59. *   longueurs des segments adequats
  60.       incnel = 2000
  61. *
  62. *  Creation de la raideur des conditions de contact
  63. *  Remplissage de l'entete commun
  64. *
  65.       nrigel = 1
  66.       segini mrigid
  67.       ichole = 0
  68.       imgeo1 = 0
  69.       imgeo2 = 0
  70.       isupeq = 0
  71.       iforig = ifour
  72.       coerig(1)=1.d0
  73.       mtymat='RIGIDITE'
  74. *
  75. *     MCHAML materiau associe au MMODEL     
  76.       MELVA1 = 0
  77.       MELVA2 = 0
  78.       IF (MCHEL1.NE.0) THEN
  79.         SEGACT, MCHEL1
  80. *  recherche rapide d'un maillage correspondant dans le mchaml
  81.         DO 210 n2 = 1,mchel1.imache(/1)
  82. *        write(ioimp,*) ' n2 imache(n2) ipt1 ',n2,mchel1.imache(n2),ipt1
  83.           if (mchel1.imache(n2).eq.ipt1) goto 220
  84.  210    continue
  85.  
  86.         goto 230
  87.  220    continue
  88.         MCHAM1 = MCHEL1.ICHAML(n2)
  89.         SEGACT, MCHAM1
  90.         NBCO = MCHAM1.NOMCHE(/2)
  91.         DO JCMP = 1,NBCO
  92.           IF (MCHAM1.NOMCHE(JCMP).EQ.'JEU') THEN
  93.             MELVA1 = MCHAM1.IELVAL(JCMP)
  94.             SEGACT, MELVA1
  95.             NELJ = MELVA1.VELCHE(/2)
  96.             NPTELJ = min(MELVA1.VELCHE(/1),4)
  97. C            WRITE(*,*) 'NELJ NPTELJ ',NELJ,NPTELJ
  98.           ENDIF
  99.           IF (MCHAM1.NOMCHE(JCMP).EQ.'ADHE') THEN
  100.             MELVA2 = MCHAM1.IELVAL(JCMP)
  101.             SEGACT, MELVA2
  102.             NELA = MELVA2.VELCHE(/2)
  103.             NPTELA = min(MELVA2.VELCHE(/1),4)
  104.           ENDIF
  105.         ENDDO
  106.       ENDIF
  107.  230  continue
  108. *
  109. *     Creation du chpoint de depi
  110. *
  111.       nat=1
  112.       nsoupo=1
  113.       segini mchpoi
  114.       mtypoi='DEPIMP'
  115.       mochde='engendré par impose'
  116.       ifopoi=ifour
  117.       jattri(1)=2
  118. *
  119.       nc=2
  120.       IF (MELVA2.NE.0) nc=3
  121.       segini msoupo
  122.       ipchp(1)=msoupo
  123.       nocomp(1)='FLX '
  124.       nocomp(2)='SCAL'
  125.       IF (MELVA2.NE.0) THEN 
  126.         nocomp(3)='FADH'
  127.       ENDIF
  128. *
  129.       nbnn  =1
  130.       nbelem=incnel
  131.       nbref =0
  132.       nbsous=0
  133.       segini ipt7
  134.       ipt7.itypel=1
  135.       igeoc=ipt7
  136. *
  137.       n=incnel
  138.       segini mpoval
  139.       ipoval=mpoval
  140. *      
  141. *  Calcul des directions moyennes
  142. **    segact mcoord*mod
  143.       segini mfopa2
  144. *      write (6,*) ' impos2 iel1 ',nbel1,idimp1
  145.       do 820 iel6=1,nbel6
  146. *
  147.         ip1=ipt6.num(1,iel6)
  148.         ip2=ipt6.num(2,iel6)
  149.         ipv1 = (ip1-1)*idimp1
  150.         xp1 = xcoor(ipv1+1)
  151.         yp1 = xcoor(ipv1+2)
  152.         ipv2 = (ip2-1)*idimp1
  153.         xp2 = xcoor(ipv2+1)
  154.         yp2 = xcoor(ipv2+2)
  155. *
  156. *    normale a la droite (12)
  157. *
  158.         x12 = xp2 - xp1
  159.         y12 = yp2 - yp1
  160.         xn =  -y12
  161.         yn =   x12
  162.         sn = sqrt (xn*xn + yn*yn)
  163.         sn = max(xpetit,sn)
  164.         xn = xn/sn
  165.         yn = yn/sn
  166.         xms(ip1)=xms(ip1)+xn
  167.         yms(ip1)=yms(ip1)+yn
  168.         xms(ip2)=xms(ip2)+xn
  169.         yms(ip2)=yms(ip2)+yn
  170.  820  continue
  171.       do 822 iel8=1,nbel8
  172. *
  173.         ip1=ipt8.num(1,iel8)
  174.         ip2=ipt8.num(2,iel8)
  175.         ipv1 = (ip1-1)*idimp1
  176.         xp1 = xcoor(ipv1+1)
  177.         yp1 = xcoor(ipv1+2)
  178.         ipv2 = (ip2-1)*idimp1
  179.         xp2 = xcoor(ipv2+1)
  180.         yp2 = xcoor(ipv2+2)
  181. *
  182. *    normale a la droite (12)
  183. *
  184.         x12 = xp2 - xp1
  185.         y12 = yp2 - yp1
  186.         xn =  -y12
  187.         yn =   x12
  188.         sn = sqrt (xn*xn + yn*yn)
  189.         sn = max(xpetit,sn)
  190.         xn = xn/sn
  191.         yn = yn/sn
  192.         xms(ip1)=xms(ip1)+xn
  193.         yms(ip1)=yms(ip1)+yn
  194.         xms(ip2)=xms(ip2)+xn
  195.         yms(ip2)=yms(ip2)+yn
  196.  822  continue
  197.         do 821 ip = 1,nbpts
  198.          sn = xms(ip)*xms(ip)+yms(ip)*yms(ip)
  199.          sn = max(sqrt(abs(sn)),xpetit*1d10)
  200.          xms(ip)=xms(ip)/sn
  201.          yms(ip)=yms(ip)/sn
  202.  821  continue
  203. *
  204. *
  205. *  Nombre de noeuds dans le chpoint (en totalite) : ipt7 et mpoval
  206.       nelch = 0
  207. *
  208. *=======================================================================
  209. * Formulation "faible" du contact : relation segment-segment (type 0)
  210. *=======================================================================
  211. * Nouvelle formulation, une seule relation par element maitre
  212. * Il faut donc reunir les relations portant sur le meme multiplicateur
  213. ************************************************************************
  214.  
  215. *   itcont 2 formulation faible
  216.       if (itcont.ne.2) goto 1500
  217. *
  218. *  Creation du meleme associe a la relation
  219.       nbnn  =5
  220.       nbelem=incnel
  221.       nbsous=0
  222.       nbref =0
  223.       segini meleme
  224.       itypel=22
  225.       irigel(1,nrigel) = meleme
  226. *
  227. *  Creation du descriptif commun a toutes les raideurs
  228. *
  229.       nligrp=9
  230.       nligrd=nligrp
  231.       segini,descr
  232.       lisinc(1)='LX  '
  233.       lisdua(1)='FLX '
  234.       noelep(1)=1
  235.       noeled(1)=1
  236.       do 100 i = 2, nligrp, 2
  237.         lisinc(i  )=modepl(imo)
  238.         lisinc(i+1)=modepl(imo+1)
  239.         lisdua(i  )=moforc(imo)
  240.         lisdua(i+1)=moforc(imo+1)
  241.         noelep(i  )=(i+2)/2
  242.         noelep(i+1)=noelep(i)
  243.         noeled(i  )=noelep(i)
  244.         noeled(i+1)=noelep(i)
  245.  100  continue
  246.       segdes,descr
  247.       irigel(3,nrigel) = descr
  248. *
  249. *     creation du segment xmatri
  250. *
  251.       nelrig=incnel
  252.       segini xmatri
  253.       irigel(4,nrigel) = xmatri
  254. *
  255. *     ce qu'on cree est unilateral
  256. *
  257.       irigel(6,nrigel) = 1
  258. *
  259. *     ce qu'on cree est symetrique
  260. *
  261.       irigel(7,nrigel) = 0
  262. *
  263. *  Nombre d'elements crees dans meleme=irigel(nrigel,1), ipt7 et mpoval
  264.       nelri0 = 0
  265. *
  266. *  boucle sur le maillage support
  267. *
  268.       do 111 iel8=1,nbel8
  269. *
  270.         xjr = 0d0
  271.         if (MELVA1.ne.0) then
  272.           nel1 = min (iel8,nelj)
  273.           xjr  = melva1.velche(nptelj,nel1)
  274.         endif
  275.         ip1 = ipt8.num(1,iel8)
  276.         ip2 = ipt8.num(2,iel8)
  277.         ipv = (ip1-1)*idimp1
  278.         xp1 = xcoor(ipv+1)
  279.         yp1 = xcoor(ipv+2)
  280.         ipv = (ip2-1)*idimp1
  281.         xp2 = xcoor(ipv+1)
  282.         yp2 = xcoor(ipv+2)
  283.         xp12 = xp2 - xp1
  284.         yp12 = yp2 - yp1
  285.         d12 = ((xp12**2)+(yp12**2))**0.5
  286. *
  287.       do 110 iel6 = 1, nbel6
  288.         ip3 = ipt6.num(1,iel6)
  289.         ip4 = ipt6.num(2,iel6)
  290. *        write(ioimp,*) iel,ip1,ip2,ip3,ip4,ipt1.num(1,iel)
  291. *  verification (provisoire) que pas de noeuds doubles
  292.         if (ip1.eq.ip3) goto 110
  293.         if (ip1.eq.ip4) goto 110
  294.         if (ip2.eq.ip3) goto 110
  295.         if (ip2.eq.ip4) goto 110
  296. *
  297.         ipv = (ip3-1)*idimp1
  298.         xp3 = xcoor(ipv+1)
  299.         yp3 = xcoor(ipv+2)
  300.         ipv = (ip4-1)*idimp1
  301.         xp4 = xcoor(ipv+1)
  302.         yp4 = xcoor(ipv+2)
  303.         xp34 = xp4 - xp3
  304.         yp34 = yp4 - yp3
  305.         d34 = ((xp34**2)+(yp34**2))**0.5
  306. *
  307. *   orientations respectives correctes des 2 segments :
  308. *   "normales de sens opposes" = produit scalaire negatif ou nul
  309.         scal = xp12*xp34 + yp12*yp34
  310.         xl12 = sqrt(xp12**2+yp12**2)
  311.         xl34 = sqrt(xp34**2+yp34**2)
  312. *       if (scal.gt.-xl12*xl34*0.5) goto 110
  313.         if (scal.gt.0.d0) goto 110
  314. *
  315. *   critere d'acceptation de l'élément :
  316. *   angles des diagonales
  317.         xl13 = sqrt((xp3-xp1)**2+(yp3-yp1)**2)
  318.         xl14 = sqrt((xp4-xp1)**2+(yp4-yp1)**2)
  319.         xl23 = sqrt((xp3-xp2)**2+(yp3-yp2)**2)
  320.         xl24 = sqrt((xp4-xp2)**2+(yp4-yp2)**2)
  321.         sca312 = (xp3-xp1)*(xp3-xp2)+(yp3-yp1)*(yp3-yp2)
  322.         sca412 = (xp4-xp1)*(xp4-xp2)+(yp4-yp1)*(yp4-yp2)
  323.         sca134 = (xp1-xp3)*(xp1-xp4)+(yp1-yp3)*(yp1-yp4)
  324.         sca234 = (xp2-xp3)*(xp2-xp4)+(yp2-yp3)*(yp2-yp4)
  325. *        write(ioimp,*) sca312/(xl13*xl23),sca412/(xl14*xl24),
  326. *    &                  sca134/(xl13*xl14),sca234/(xl23*xl24)
  327. **      if (sca312/(xl13*xl23).gt.0.50.and.
  328. **   >      sca412/(xl14*xl24).gt.0.50.and.
  329. **   >      sca134/(xl13*xl14).gt.0.50.and.
  330. **   >      sca234/(xl23*xl24).gt.0.50)  goto 110
  331.  
  332. *   nouveau critere acceptation
  333. *   pts dans un cercle centree sur 1-2 aggrandi
  334.         xp1e=xp1-(xp2-xp1)
  335.         yp1e=yp1-(yp2-yp1)
  336.         xp2e=xp2+(xp2-xp1)
  337.         yp2e=yp2+(yp2-yp1)
  338. *
  339. *       Tenir compte du jeu dans le critere de selection
  340.         xpp3=xp3
  341.         ypp3=yp3
  342.         xpp4=xp4
  343.         ypp4=yp4
  344.         if (MELVA1.ne.0) then
  345.           xnorm=max(xl12,xpetit)
  346.           xn =  yp12/xnorm
  347.           yn = -xp12/xnorm
  348.           xjrxn = xn * xjr
  349.           xjryn = yn * xjr
  350.           xpp3=xpp3-xjrxn
  351.           ypp3=ypp3-xjryn
  352.           xpp4=xpp4-xjrxn
  353.           ypp4=ypp4-xjryn
  354.         endif
  355. *
  356.         sca312 = (xpp3-xp1e)*(xpp3-xp2e)+(ypp3-yp1e)*(ypp3-yp2e)
  357.         sca412 = (xpp4-xp1e)*(xpp4-xp2e)+(ypp4-yp1e)*(ypp4-yp2e)
  358.         if (sca312.gt.0.d0.and.
  359.      >      sca412.gt.0.d0) goto 110
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369. *
  370. * Quelle est la relation ???
  371. *
  372. *   direction du contact unitaire
  373. *
  374.         xr = xp12 - xp34
  375.         yr = yp12 - yp34
  376.         sr = sqrt(xr*xr + yr*yr)
  377.         xr = xr/sr
  378.         yr = yr/sr
  379. *        write(ioimp,*) 'direction contact',xr,yr,yr,-xr
  380. *
  381. *   projection des points sur la direction du contact
  382. *
  383.         xl1 = xp1*xr + yp1*yr
  384.         xl2 = xp2*xr + yp2*yr
  385.         xl3 = xp3*xr + yp3*yr
  386.         xl4 = xp4*xr + yp4*yr
  387.  
  388. *        write(ioimp,*) 'projection pts sur contact',xl1,xl2,xl3,xl4
  389. *
  390. *   calcul de l'intersection des projections
  391. *
  392.         xm1 = min(xl1,xl2)
  393.         xm2 = max(xl1,xl2)
  394.         xm3 = min(xl3,xl4)
  395.         xm4 = max(xl3,xl4)
  396. *        write(ioimp,*) ' xmi',xm1,xm2,xm3,xm4
  397. *
  398. *  critere de precision sur l'intersection
  399. *
  400.         xcr = min(xm2-xm1,xm4-xm3)*(1.d-10)
  401. *
  402. *  taille de l'intersection
  403. *
  404.         xi = max(xm1,xm3)
  405.         xj = min(xm2,xm4)
  406.         xl = xj - xi
  407. *       write(ioimp,*) ' intersection',xi,xj,xl,xcr
  408. *
  409. *              write (6,*) ' impos2 ',ip1,ip2,ip3,ip4,xcr,xl
  410. *       if (xl.le.xcr) goto 110
  411.         if (xl.le.0.d0) goto 110
  412. *
  413. *   distance des points a leur projection
  414. *
  415.         d1 = (xp1-xl1*xr)*yr-(yp1-xl1*yr)*xr
  416.         d2 = (xp2-xl2*xr)*yr-(yp2-xl2*yr)*xr
  417.         d3 = (xp3-xl3*xr)*yr-(yp3-xl3*yr)*xr
  418.         d4 = (xp4-xl4*xr)*yr-(yp4-xl4*yr)*xr
  419. *
  420. *  coordonnées paramétriques de l'intersection sur les segments 1-2 et 3-4
  421. *
  422.         xi2 = (xi-xl1) / (xl2-xl1)
  423.         xi1 = (xl2-xi) / (xl2-xl1)
  424.  
  425.         xj2 = (xj-xl1) / (xl2-xl1)
  426.         xj1 = (xl2-xj) / (xl2-xl1)
  427.  
  428.         xi4 = (xi-xl3) / (xl3-xl4)
  429.         xi3 = (xl4-xi) / (xl3-xl4)
  430.  
  431.         xj4 = (xj-xl3) / (xl3-xl4)
  432.         xj3 = (xl4-xj) / (xl3-xl4)
  433. *
  434. *        write(ioimp,*) ' xi1 xi2 xi3 xi4 xj1 xj2 xj3 xj4 xl '
  435. *        write(ioimp,*)  xi1,xi2,xi3,xi4,xj1,xj2,xj3,xj4,xl
  436. *        write(ioimp,*) ' d1 d2 d3 d4 ',d1,d2,d3,d4
  437. *
  438. * surface actuelle
  439. *
  440.         sc = ((xi1+xj1)*d1+(xi2+xj2)*d2+(xi3+xj3)*d3+(xi4+xj4)*d4)*0.5
  441. *        write(ioimp,*) 'Surface actuelle :',sc
  442. *
  443. * on a un element ou imposer la relation a ajouter
  444. *
  445.         nelri0 = nelri0 + 1
  446.         nelch  = nelch +1
  447. *
  448. *   on ajuste les differents segments
  449. *
  450.         if (nelri0.gt.nelrig) then
  451.           nelrig = nelrig + incnel
  452.           segadj,xmatri
  453.           nbelem = nbelem + incnel
  454.           nbnn = 5
  455.           segadj,meleme
  456.           nbnn = 1
  457.           segadj,ipt7
  458.           n = n + incnel
  459.           segadj,mpoval
  460.         endif
  461. *
  462. *       Choix du mult de Lagrange qui porte la condition
  463. *       -> l'element le plus grand
  464.         imult=ipt1.num(1,iel8)
  465.         ielt=iel8
  466.         if (d12.lt.d34) then
  467.           imult=ipt1.num(1,nbel8+iel6)
  468.           ielt=nbel8+iel6
  469.         endif
  470. *
  471. *   on range dans le meleme
  472. *
  473.         num(1,nelri0) = imult
  474.         num(2,nelri0) = ip1
  475.         num(3,nelri0) = ip2
  476.         num(4,nelri0) = ip3
  477.         num(5,nelri0) = ip4
  478.         icolor(nelri0)=1
  479. *
  480. *   on remplit le xmatri
  481. *
  482. *  lambda
  483.         re(1,1,nelri0)= 0.d0
  484. *  ip1
  485.         re(2,1,nelri0)=  yr * (xi1+xj1) * 0.5  * xl
  486.         re(3,1,nelri0)= -xr * (xi1+xj1) * 0.5  * xl
  487. *  ip2
  488.         re(4,1,nelri0)=  yr * (xi2+xj2) * 0.5  * xl
  489.         re(5,1,nelri0)= -xr * (xi2+xj2) * 0.5  * xl
  490. *  ip3
  491.         re(6,1,nelri0)=  yr * (xi3+xj3) * 0.5  * xl
  492.         re(7,1,nelri0)= -xr * (xi3+xj3) * 0.5  * xl
  493. *  ip4
  494.         re(8,1,nelri0)=  yr * (xi4+xj4) * 0.5  * xl
  495.         re(9,1,nelri0)= -xr * (xi4+xj4) * 0.5  * xl
  496. *
  497. *  on transpose
  498.         do 120 ic = 2, nligrp
  499.           re(1,ic,nelri0)=re(ic,1,nelri0)
  500.  120    continue
  501. *
  502. *  Le reste est nul
  503. *
  504. *  remplissage du champoint de depi et du maillage support
  505. *
  506.         ipt7.num(1,nelch)=imult
  507.         vpocha(nelch,1) = (-sc - xjr) * xl
  508.         vpocha(nelch,2) = xl
  509.         IF (MELVA2.NE.0) THEN  
  510.           NEL2 = min (ielt,NELA)
  511.           VPOCHA(nelch,3) = max(MELVA2.VELCHE(NPTELA,NEL2),0.D0)*xl
  512.         ENDIF  
  513. *     write (6,*) 'impos2 vpocha re ',vpocha(nelch,1),
  514. *    >  (re(ip,1,nelri0),ip=2,9)
  515. *
  516.  110   continue
  517.  111   continue
  518.  
  519. *     write (ioimp,*) ' nb relation type 0 ',nbelem,n,nelri0
  520.  
  521. * Ajustement au plus juste puis desactivation des segments lies
  522. * a la rigidite du type 0
  523.       if (nelri0.ne.nelrig) then
  524.         nelrig = nelri0
  525.         segadj,xmatri
  526.         nbelem = nelri0
  527.         nbnn = 5
  528.         segadj,meleme
  529.       endif
  530.       segdes,xmatri
  531. *
  532. * S'il n'y a pas d'elements en contact, alors pas de relation unilaterale
  533. *     if (nelri0.eq.0) irigel(6,nrigel)=0
  534.  
  535.  
  536.  
  537.  
  538.  
  539.  
  540.       GOTO 1000
  541.  
  542. *
  543. *=======================================================================
  544. * Formulation "mortar" du contact 
  545. *=======================================================================
  546. * 
  547. * 
  548. ************************************************************************
  549.  
  550.  1500 continue
  551.  
  552. *   itcont 3 formulation mortar
  553.       if (itcont.ne.4) goto 500
  554. *
  555.  
  556. *  Creation du meleme associe a la relation
  557.       nbnn  =5
  558.       nbelem=incnel
  559.       nbsous=0
  560.       nbref =0
  561.       segini meleme
  562.       itypel=22
  563.       irigel(1,nrigel) = meleme
  564. *
  565. *  Creation du descriptif commun a toutes les raideurs
  566. *
  567.       nligrp=9
  568.       nligrd=nligrp
  569.       segini,descr
  570.       lisinc(1)='LX  '
  571.       lisdua(1)='FLX '
  572.       noelep(1)=1
  573.       noeled(1)=1
  574.       do 300 i = 2, nligrp, 2
  575.         lisinc(i  )=modepl(imo)
  576.         lisinc(i+1)=modepl(imo+1)
  577.         lisdua(i  )=moforc(imo)
  578.         lisdua(i+1)=moforc(imo+1)
  579.         noelep(i  )=(i+2)/2
  580.         noelep(i+1)=noelep(i)
  581.         noeled(i  )=noelep(i)
  582.         noeled(i+1)=noelep(i)
  583.  300  continue
  584.       segdes,descr
  585.       irigel(3,nrigel) = descr
  586. *
  587. *     creation du segment xmatri
  588. *
  589.       nelrig=incnel
  590.       segini xmatri
  591.       irigel(4,nrigel) = xmatri
  592. *
  593. *     ce qu'on cree est unilateral
  594. *
  595.       irigel(6,nrigel) = 1
  596. *
  597. *     ce qu'on cree est symetrique
  598. *
  599.       irigel(7,nrigel) = 0
  600. *
  601. *  Nombre d'elements crees dans meleme=irigel(nrigel,1), ipt7 et mpoval
  602.       nelri0 = 0
  603. *
  604. *  boucle sur le maillage non mortar
  605. *
  606.       do 311 iel8 = 1, nbel8
  607. *
  608. *  Si jeu :
  609.         xjr = 0d0
  610.         if (MELVA1.ne.0) then
  611.            nel8 = min (iel8,nelj)
  612.            xjr = melva1.velche(nptelj,nel8)
  613.         endif
  614.         ip1 = ipt8.num(1,iel8)
  615.         ip2 = ipt8.num(2,iel8)
  616.         ipv = (ip1-1)*idimp1
  617.         xp1 = xcoor(ipv+1)
  618.         yp1 = xcoor(ipv+2)
  619.         ipv = (ip2-1)*idimp1
  620.         xp2 = xcoor(ipv+1)
  621.         yp2 = xcoor(ipv+2)
  622.         xp12 = xp2 - xp1
  623.         yp12 = yp2 - yp1
  624. *
  625.       do 310 iel6 = 1, nbel6
  626.         ip3 = ipt6.num(1,iel6)
  627.         ip4 = ipt6.num(2,iel6)
  628. *        write(ioimp,*) iel,ip1,ip2,ip3,ip4,ipt1.num(1,iel)
  629. *  verification (provisoire) que pas de noeuds doubles
  630.         if (ip1.eq.ip3) goto 310
  631.         if (ip1.eq.ip4) goto 310
  632.         if (ip2.eq.ip3) goto 310
  633.         if (ip2.eq.ip4) goto 310
  634. *
  635.         ipv = (ip3-1)*idimp1
  636.         xp3 = xcoor(ipv+1)
  637.         yp3 = xcoor(ipv+2)
  638.         ipv = (ip4-1)*idimp1
  639.         xp4 = xcoor(ipv+1)
  640.         yp4 = xcoor(ipv+2)
  641.         xp34 = xp4 - xp3
  642.         yp34 = yp4 - yp3
  643. *
  644. *   orientations respectives correctes des 2 segments :
  645. *   "normales de sens opposes" = produit scalaire negatif ou nul
  646.         scal = xp12*xp34 + yp12*yp34
  647.         xl12 = sqrt(xp12**2+yp12**2)
  648.         xl34 = sqrt(xp34**2+yp34**2)
  649.         if (scal.gt.0.d0) goto 310
  650. *
  651. *   nouveau critere acceptation
  652. *   pts dans un cercle centree sur 1-2 aggrandi
  653.         xp1e=xp1-(xp2-xp1)
  654.         yp1e=yp1-(yp2-yp1)
  655.         xp2e=xp2+(xp2-xp1)
  656.         yp2e=yp2+(yp2-yp1)
  657. *
  658. *       Tenir compte du jeu dans le critere de selection
  659.         xpp3=xp3
  660.         ypp3=yp3
  661.         xpp4=xp4
  662.         ypp4=yp4
  663. *   Si jeu :
  664.         if (MELVA1.ne.0) then
  665.           xnorm=max(xl12,xpetit)
  666.           xn =  yp12/xnorm
  667.           yn = -xp12/xnorm
  668.           xjrxn = xn * xjr
  669.           xjryn = yn * xjr
  670.           xpp3=xpp3-xjrxn
  671.           ypp3=ypp3-xjryn
  672.           xpp4=xpp4-xjrxn
  673.           ypp4=ypp4-xjryn
  674.         endif
  675. *
  676.         sca312 = (xpp3-xp1e)*(xpp3-xp2e)+(ypp3-yp1e)*(ypp3-yp2e)
  677.         sca412 = (xpp4-xp1e)*(xpp4-xp2e)+(ypp4-yp1e)*(ypp4-yp2e)
  678.         if (sca312.gt.0.d0.and.
  679.      >      sca412.gt.0.d0) goto 310
  680.  
  681. *   Vecteurs unitaires tangent et normal (element Mortar)
  682.         xt34 = xp34/xl34
  683.         yt34 = yp34/xl34
  684.         xn34 = yt34
  685.         yn34 = -xt34
  686. *
  687. *   Projection des points non-Mortar sur le segment Mortar
  688. *
  689.         xl1 = (xp1-xp3)*xt34 + (yp1-yp3)*yt34
  690.         xl2 = (xp2-xp3)*xt34 + (yp2-yp3)*yt34
  691. *
  692. *   Intersection des projections et taille de l'intersection
  693. *
  694.         xm1 = min(xl1,xl2)
  695.         xm2 = max(xl1,xl2)
  696.         xi  = max(xm1,0.D0)
  697.         xj  = min(xm2,xl34)
  698.         xlmo  = xj - xi
  699. *       if (xlmo.le.0.d0) goto 310
  700.         if (xlmo.le.1.D-10) goto 310
  701. *
  702. *   Coordonnees parametriques de l'intersection des segments
  703. *
  704. *   - sur le segment 1-2 
  705.         xtest12 = xm2 - xm1
  706.         xi1 = (xi-xl2) / xtest12
  707.         xi2 = 1.D0 - xi1
  708.  
  709.         xj1 = (xj-xl2) / xtest12
  710.         xj2 = 1.D0 - xj1
  711. *
  712. *   - sur le segment 3-4 
  713.         xi4 = xi / xl34
  714.         xi3 = 1.D0 - xi4
  715.  
  716.         xj4 = xj / xl34
  717.         xj3 = 1.D0 - xj4
  718. *
  719. *   Coordonnees des points limites de la surface de contact
  720. *
  721. *   - sur le segment 1-2 
  722.         xnm1 = xp1 + xi2*xp12
  723.         ynm1 = yp1 + xi2*yp12
  724.         xnm2 = xp1 + xj2*xp12
  725.         ynm2 = yp1 + xj2*yp12
  726.         xnm12 = xnm2 - xnm1
  727.         ynm12 = ynm2 - ynm1
  728.         xlnm  = sqrt(xnm12**2 + ynm12**2) 
  729. *
  730. *   - sur le segment 3-4 
  731.         xmo1 = xp3 + xi4*xp34
  732.         ymo1 = yp3 + xi4*yp34
  733.         xmo2 = xp3 + xj4*xp34
  734.         ymo2 = yp3 + xj4*yp34
  735.         xmo12 = xmo2 - xmo1
  736.         ymo12 = ymo2 - ymo1
  737. *
  738. *   Points de Gauss a positionner sur l'intersection cote Non Mortar
  739. *
  740.         xref1 = 0.5 - 0.5*(1.d0/sqrt(3.d0))
  741.         xref2 = 0.5 + 0.5*(1.d0/sqrt(3.d0))
  742. *
  743. *   Coordonnees des points de Gauss sur Non Mortar
  744. *
  745.         xksi1 = xnm1+xref1*xnm12
  746.         yksi1 = ynm1+xref1*ynm12
  747.         xksi2 = xnm1+xref2*xnm12
  748.         yksi2 = ynm1+xref2*ynm12
  749. C
  750.         xnmg1 = (xksi1-xnm1)*xt34 + (yksi1-ynm1)*yt34
  751.         xnmg2 = (xksi2-xnm1)*xt34 + (yksi2-ynm1)*yt34
  752.         zksi1 = xnmg1 / xlmo
  753.         zksi2 = xnmg2 / xlmo
  754. *
  755. *   Projection des points de Gauss Non Mortar sur Mortar
  756. *
  757.         xmog1 = (xksi1-xmo1)*xt34 + (yksi1-ymo1)*yt34
  758.         xmog2 = (xksi2-xmo1)*xt34 + (yksi2-ymo1)*yt34
  759.         zeta1 = xmog1 / xlmo
  760.         zeta2 = xmog2 / xlmo
  761. C
  762.         xtau1 = xmo1+zeta1*xmo12
  763.         ytau1 = ymo1+zeta1*ymo12
  764.         xtau2 = xmo1+zeta2*xmo12
  765.         ytau2 = ymo1+zeta2*ymo12
  766. C
  767.         xiG(1)   = ((xi1*xl12) + (zksi1*xlnm)) / xl12
  768.         xiG(2)   = ((xi1*xl12) + (zksi2*xlnm)) / xl12
  769.         zetaG(1) = ((xi4*xl34) + (zeta1*xlmo)) / xl34
  770.         zetaG(2) = ((xi4*xl34) + (zeta2*xlmo)) / xl34
  771. C
  772.         wG(1) = 1.D0
  773.         wG(2) = 1.D0
  774. C
  775.         if (iimpi.eq.2505) then
  776.           write(*,*) 'NOUVELLE RELATION DE CONTACT ENTRE LES NOEUDS :' 
  777.           write(*,*) '==============================================='
  778.           write(*,*) '(ELEMENT ',iel8, 'DU MAILLAGE ',ipt8,')'
  779.           write(*,*) '- IP1=',ip1,'x=',xp1,'y=',yp1
  780.           write(*,*) '- IP2=',ip2,'x=',xp2,'y=',yp2
  781.           write(*,*) 'INTERSECTION DE CONTACT SUR IP1-IP2, L=',xlnm
  782.           write(*,*) '- POINT1(x,y) PT1=',xnm1,ynm1,';'
  783.           write(*,*) '- POINT2(x,y) PT2=',xnm2,ynm2,';'
  784.           write(*,*) 'COORDONNES DES POINTS DE GAUSS SUR IP1-IP2'
  785.           write(*,*) '- PDG1(x,y) PG1=',xksi1,yksi1,';'
  786.           write(*,*) '- PDG2(x,y) PG2=',xksi2,yksi2,';'
  787.           write(*,*) 'COEF PDG1 IP1=',xiG(1),' IP2=',(1.D0-xiG(1))
  788.           write(*,*) 'COEF PDG2 IP1=',xiG(2),' IP2=',(1.D0-xiG(2))
  789. C
  790.           write(*,*) '(ELEMENT ',iel6, 'DU MAILLAGE ',ipt6,')'
  791.           write(*,*) '- IP3=',ip3,'x=',xp3,'y=',yp3
  792.           write(*,*) '- IP4=',ip4,'x=',xp4,'y=',yp4
  793.           write(*,*) '- VECTEUR NORMAL ', xn34,yn34
  794.           write(*,*) 'INTERSECTION DE CONTACT SUR IP3-IP4, L=',xlmo
  795.           write(*,*) '- POINT1(x,y) PT1=',xmo1,ymo1,';'
  796.           write(*,*) '- POINT1(x,y) PT2=',xmo2,ymo2,';'
  797.           write(*,*) 'COORDONNES DES POINTS DE GAUSS SUR IP1-IP2'
  798.           write(*,*) '- PDG1(x,y) PG1=',xtau1,ytau1,';'
  799.           write(*,*) '- PDG2(x,y) PG2=',xtau2,ytau2,';'
  800.           write(*,*) 'COEF PDG1 IP3=',(1.D0-zetaG(1)),' IP4=',zetaG(1)
  801.           write(*,*) 'COEF PDG2 IP3=',(1.D0-zetaG(2)),' IP4=',zetaG(2)
  802.         endif
  803. *
  804. * on a un element ou imposer la relation a ajouter
  805. *
  806.         nelri0 = nelri0 + 2
  807.         nelch  = nelch  + 2
  808. *
  809. *   on ajuste les differents segments
  810. *
  811.         if (nelri0.gt.nelrig) then
  812.           nelrig = nelrig + incnel
  813.           segadj,xmatri
  814.           nbelem = nbelem + incnel
  815.           nbnn = 5
  816.           segadj,meleme
  817.           nbnn = 1
  818.           segadj,ipt7
  819.           n = n + incnel
  820.           segadj,mpoval
  821.         endif
  822. *
  823. *   on range dans le meleme
  824. *
  825. *   Il faut retrouver dans ipt1 les noeuds support des mult de Lag.
  826.         ilamb1 = 0
  827.         ilamb2 = 0
  828.         do 330 iel1 = 1, nbel1
  829.           if ((ipt1.num(2,iel1)).eq.ip1) then
  830.             ilamb1 = ipt1.num(1,iel1)
  831.           endif
  832.           if ((ipt1.num(2,iel1)).eq.ip2) then
  833.             ilamb2 = ipt1.num(1,iel1)
  834.           endif
  835.           if (ilamb1.ne.0.and.ilamb2.ne.0) GOTO 331
  836.  330    continue
  837.  331    continue
  838. CCCC    WRITE(*,*) 'RELATION ',ip1,ip2,ilamb1,ilamb2,ip3,ip4
  839. *
  840.         num(1,nelri0-1) = ilamb1
  841.         num(2,nelri0-1) = ip1
  842.         num(3,nelri0-1) = ip2
  843.         num(4,nelri0-1) = ip3
  844.         num(5,nelri0-1) = ip4
  845. *
  846.         num(1,nelri0) = ilamb2
  847.         num(2,nelri0) = ip1
  848.         num(3,nelri0) = ip2
  849.         num(4,nelri0) = ip3
  850.         num(5,nelri0) = ip4
  851. *
  852. **      icolor(nelri0)=2
  853.  
  854. *   on initialise le xmatri
  855. *
  856. *  Matrice elementaire associee au multiplicateur lambda1
  857.         re(1,1,nelri0-1)= 0.d0
  858. *  ip1 - lambda1
  859.         re(2,1,nelri0-1)= 0.d0
  860.         re(3,1,nelri0-1)= 0.d0
  861. *  ip2 - lambda1
  862.         re(4,1,nelri0-1)= 0.d0
  863.         re(5,1,nelri0-1)= 0.d0 
  864. *  ip3 - lambda1
  865.         re(6,1,nelri0-1)= 0.d0
  866.         re(7,1,nelri0-1)= 0.d0
  867. *  ip4 - lambda1
  868.         re(8,1,nelri0-1)= 0.d0 
  869.         re(9,1,nelri0-1)= 0.d0
  870. *
  871. *  Matrice elementaire associee au multiplicateur lambda2
  872.         re(1,1,nelri0)= 0.d0
  873. *  ip1 - lambda2
  874.         re(2,1,nelri0)= 0.d0 
  875.         re(3,1,nelri0)= 0.d0
  876. *  ip2 - lambda2
  877.         re(4,1,nelri0)= 0.d0 
  878.         re(5,1,nelri0)= 0.d0
  879. *  ip3 - lambda2
  880.         re(6,1,nelri0)= 0.d0 
  881.         re(7,1,nelri0)= 0.d0
  882. *  ip4 - lambda2
  883.         re(8,1,nelri0)= 0.d0 
  884.         re(9,1,nelri0)= 0.d0 
  885. *
  886. *  CHPOINT depi
  887.         vpocha(nelch-1,1) = 0.d0
  888.         vpocha(nelch,1)   = 0.d0
  889. * 
  890.         do 350 iGauss = 1,2
  891. *   On recupere les valeurs des pts de Gauss et des poids
  892.           xiGi = xiG(iGauss)
  893.           zetaGi = zetaG(iGauss)
  894.           wGi = wG(iGauss)
  895. *
  896. *   On evalue les fonctions d'interpolation au pt de Gauss
  897.           fM1Gi = xiGi
  898.           fM2Gi = 1.D0 - xiGi
  899. *
  900.           fN1Gin = xiGi
  901.           fN2Gin = 1.D0 - xiGi
  902.           fN1Gim = 1.D0 - zetaGi
  903.           fN2Gim = zetaGi
  904. *
  905. *   on remplit le xmatri
  906. *
  907. *  Matrice elementaire associee au multiplicateur lambda1
  908. *         re(1,1,nelri0-1)= 0.d0
  909. *  ip1 - lambda1
  910.           re(2,1,nelri0-1)=re(2,1,nelri0-1) - wGi*fM1Gi*xn34*fN1Gin*xlnm
  911.           re(3,1,nelri0-1)=re(3,1,nelri0-1) - wGi*fM1Gi*yn34*fN1Gin*xlnm
  912. *  ip2 - lambda1
  913.           re(4,1,nelri0-1)=re(4,1,nelri0-1) - wGi*fM1Gi*xn34*fN2Gin*xlnm
  914.           re(5,1,nelri0-1)=re(5,1,nelri0-1) - wGi*fM1Gi*yn34*fN2Gin*xlnm
  915. *  ip3 - lambda1
  916.           re(6,1,nelri0-1)=re(6,1,nelri0-1) + wGi*fM1Gi*xn34*fN1Gim*xlnm
  917.           re(7,1,nelri0-1)=re(7,1,nelri0-1) + wGi*fM1Gi*yn34*fN1Gim*xlnm
  918. *  ip4 - lambda1
  919.           re(8,1,nelri0-1)=re(8,1,nelri0-1) + wGi*fM1Gi*xn34*fN2Gim*xlnm
  920.           re(9,1,nelri0-1)=re(9,1,nelri0-1) + wGi*fM1Gi*yn34*fN2Gim*xlnm
  921. *
  922. *  Matrice elementaire associee au multiplicateur lambda2
  923. *         re(1,2,nelri0)= 0.d0
  924. *  ip1 - lambda2
  925.           re(2,1,nelri0)= re(2,1,nelri0) - wGi*fM2Gi*xn34*fN1Gin*xlnm
  926.           re(3,1,nelri0)= re(3,1,nelri0) - wGi*fM2Gi*yn34*fN1Gin*xlnm
  927. *  ip2 - lambda2
  928.           re(4,1,nelri0)= re(4,1,nelri0) - wGi*fM2Gi*xn34*fN2Gin*xlnm
  929.           re(5,1,nelri0)= re(5,1,nelri0) - wGi*fM2Gi*yn34*fN2Gin*xlnm
  930. *  ip3 - lambda2
  931.           re(6,1,nelri0)= re(6,1,nelri0) + wGi*fM2Gi*xn34*fN1Gim*xlnm
  932.           re(7,1,nelri0)= re(7,1,nelri0) + wGi*fM2Gi*yn34*fN1Gim*xlnm
  933. *  ip4 - lambda2
  934.           re(8,1,nelri0)= re(8,1,nelri0) + wGi*fM2Gi*xn34*fN2Gim*xlnm
  935.           re(9,1,nelri0)= re(9,1,nelri0) + wGi*fM2Gi*yn34*fN2Gim*xlnm
  936.  
  937. *  CHPOINT de jeu pour lambda1 et lamda2
  938. *
  939.           xjeu = fN1Gin*xp1+fN2Gin*xp2-(fN1Gim*xp3+fN2Gim*xp4)
  940.           yjeu = fN1Gin*yp1+fN2Gin*yp2-(fN1Gim*yp3+fN2Gim*yp4)
  941.           zjeu = xjeu*xn34 + yjeu*yn34
  942. *
  943.           ipt7.num(1,nelch-1) = ilamb1
  944.           vpocha(nelch-1,1) = vpocha(nelch-1,1) + wGi*fM1Gi*zjeu*xlnm
  945. *
  946.           ipt7.num(1,nelch) = ilamb2
  947.           vpocha(nelch,1) = vpocha(nelch,1) + wGi*fM2Gi*zjeu*xlnm
  948.  
  949.  350    continue
  950. *
  951.         vpocha(nelch-1,2) = xlnm
  952.         vpocha(nelch  ,2) = xlnm
  953. *
  954. *  on transpose
  955.         do 320 ic = 2, nligrp
  956.           re(1,ic,nelri0-1)= re(ic,1,nelri0-1)
  957.           re(1,ic,nelri0)  = re(ic,1,nelri0)
  958.  320    continue
  959. *
  960.  310   continue
  961.  311   continue
  962.  
  963. *     write (ioimp,*) ' nb relation type 0 ',nbelem,n,nelri0
  964.  
  965. * Ajustement au plus juste puis desactivation des segments lies
  966. * a la rigidite du type 0
  967.       if (nelri0.ne.nelrig) then
  968.         nelrig = nelri0
  969.         segadj,xmatri
  970.         nbelem = nelri0
  971.         nbnn = 5
  972.         segadj,meleme
  973.       endif
  974.       segdes,xmatri
  975. *
  976. * S'il n'y a pas d'elements en contact, alors pas de relation unilaterale
  977. *     if (nelri0.eq.0) irigel(6,nrigel)=0
  978.       GOTO 1000
  979.  
  980. *=======================================================================
  981. * Formulation "forte" (standard) du contact :
  982. *=======================================================================
  983.  500  continue
  984. *
  985. * Relation type 2 : noeud-segment
  986. *---------------------------------
  987. *     creation du meleme associe a la relation
  988. *
  989.       nbnn   = 4
  990.       nbelem = incnel
  991.       nbsous = 0
  992.       nbref  = 0
  993.       segini meleme
  994.       itypel=22
  995.       irigel(1,nrigel)=meleme
  996. *
  997. *     creation du descriptif commun a toutes les raideurs
  998. *
  999.       nligrp=7
  1000.       nligrd=nligrp
  1001.       segini descr
  1002.       lisinc(1)='LX  '
  1003.       lisdua(1)='FLX '
  1004.       noelep(1)=1
  1005.       noeled(1)=1
  1006.       do 510 i=2,nligrp,2
  1007.         lisinc(i  )=modepl(imo)
  1008.         lisinc(i+1)=modepl(imo+1)
  1009.         lisdua(i  )=moforc(imo)
  1010.         lisdua(i+1)=moforc(imo+1)
  1011.         noelep(i  )=(i+2)/2
  1012.         noelep(i+1)=noelep(i)
  1013.         noeled(i  )=noelep(i)
  1014.         noeled(i+1)=noelep(i)
  1015.  510  continue
  1016.       segdes,descr
  1017.       irigel(3,nrigel)=descr
  1018. *
  1019. *     creation du segment xmatri
  1020. *
  1021.       nelrig=incnel
  1022.       segini xmatri
  1023.       irigel(4,nrigel)=xmatri
  1024. *
  1025. *     ce qu'on cree est unilateral
  1026. *
  1027.       irigel(6,nrigel)=1
  1028. *
  1029. *     ce qu'on cree est symetrique
  1030. *
  1031.       irigel(7,nrigel)=0
  1032. *
  1033. *  Nombre d'elements dans la rigidite de type 2
  1034.       nelri2=0
  1035. *
  1036. *  boucle sur le maillage support
  1037. *
  1038. *     write(6,*) 'nbel6 nbel1',nbel6,nbel1
  1039.  
  1040.  
  1041.       do 519 iel6 = 1,nbel6
  1042.         ip1 = ipt6.num(1,iel6)
  1043.         ip2 = ipt6.num(2,iel6)
  1044.         ipv = (ip1-1)*idimp1
  1045.         xp1 = xcoor(ipv+1)
  1046.         yp1 = xcoor(ipv+2)
  1047.         ipv = (ip2-1)*idimp1
  1048.         xp2 = xcoor(ipv+1)
  1049.         yp2 = xcoor(ipv+2)
  1050.         x12 = xp2 - xp1
  1051.         y12 = yp2 - yp1
  1052.         sr2=x12**2+y12**2
  1053.         sr= sqrt (max(xpetit,sr2))
  1054.  
  1055.         do 520 iel1=1,nbel1
  1056.         xjr = 0.d0
  1057.         if (MELVA1.ne.0) then
  1058.           nel1 = min (iel1,nelj)
  1059.           xjr = melva1.velche(nptelj,nel1)
  1060. *        write(ioimp,*) iel,xjr,mchel1
  1061.         endif
  1062.         jp  = ipt1.num(2,iel1)
  1063. *      write(ioimp,*) iel1,ip1,ip2,jp
  1064. *
  1065. *  verification que pas relation du point sur lui meme
  1066.         if (jp.eq.ip1) goto 520
  1067.         if (jp.eq.ip2) goto 520
  1068.  
  1069.         ipv = (jp-1)*idimp1
  1070.         xp  = xcoor(ipv+1)
  1071.         yp  = xcoor(ipv+2)
  1072.  
  1073.         x1p = xp - xp1
  1074.         y1p = yp - yp1
  1075.         x2p = xp - xp2
  1076.         y2p = yp - yp2
  1077. *  distance signee de p a la ligne 1-2
  1078.         dp12s = x1p * y12 - y1p * x12
  1079. *        write(ioimp,*) 'dist. signee',dp12s
  1080. *   verif que le point est du bon cote du segment (a une tolerance de ratrapage pres)
  1081. *       if (dp12s.lt.-sr2) goto 520
  1082.  
  1083. *  verification si autre point dans le cercle de selection (tres legerement agrandi)      
  1084.         tx12=x12/sr
  1085.         ty12=y12/sr
  1086. *
  1087.         x1e = xp1-tx12*xszpre 
  1088.         y1e = yp1-ty12*xszpre 
  1089.         x2e = xp2+tx12*xszpre 
  1090.         y2e = yp2+ty12*xszpre 
  1091. *
  1092. *       Tenir compte du jeu dans le critere de selection
  1093.         xpp=xp
  1094.         ypp=yp
  1095.         if (MELVA1.ne.0) then
  1096.           xn =  ty12
  1097.           yn = -tx12
  1098.           xjrxn = xn * xjr
  1099.           xjryn = yn * xjr
  1100.           xpp=xpp-xjrxn
  1101.           ypp=ypp-xjryn
  1102.         endif
  1103.         x1pe=xpp-x1e
  1104.         y1pe=ypp-y1e
  1105.         x2pe=xpp-x2e
  1106.         y2pe=ypp-y2e
  1107. *
  1108.         d1pe = x1pe*x1pe + y1pe*y1pe
  1109.         d2pe = x2pe*x2pe + y2pe*y2pe
  1110.         if (abs(d1pe).lt.XPETIT) d1pe=1
  1111.         if (abs(d2pe).lt.XPETIT) d2pe=1
  1112.         scal = (x1pe*x2pe + y1pe*y2pe) / sqrt(d1pe*d2pe)
  1113. *       write(ioimp,*) 'cos(angle_1p2)',scal,d1pe,d2pe
  1114.         if (scal.gt.0.8) goto 520
  1115. *
  1116. * on a un point ou imposer la relation
  1117. * Quelle est la relation ???
  1118. *
  1119. *   direction de la relation
  1120. *
  1121. * initialisation avec la direction normale, calcul du point projete, puis
  1122. * iteration en reestimant la normale a partir du point projete
  1123. *
  1124.          xn =  -y12
  1125.          yn =   x12
  1126. *        sn = sqrt (max(xpetit,xn*xn + yn*yn))
  1127.          xn = xn/sr
  1128.          yn = yn/sr
  1129.  
  1130.         do iter=1,16
  1131. *   calcul du pt a mettre en relation avec ip : alpha ip2 + (1-alpha) ip1
  1132. *  projection suivant la normale
  1133.            beta=(x1p*y12-y1p*x12)
  1134.            beta=beta/((xn*y12-yn*x12))
  1135.            xpr=xp-beta*xn
  1136.            ypr=yp-beta*yn
  1137.            alpha=((xpr-xp1)*x12+(ypr-yp1)*y12)/sr2
  1138.            alpha=min(max(0.d0,alpha),1.d0)
  1139. *
  1140. * nouvelle normale normalisee
  1141.         xn=(1.-alpha)*xms(ip1)+alpha*xms(ip2)
  1142.         yn=(1.-alpha)*yms(ip1)+alpha*yms(ip2)
  1143.            sn = sqrt (max(xpetit,xn*xn + yn*yn))
  1144.          xn = xn/sn
  1145.          yn = yn/sn
  1146.         enddo
  1147. * verif dans le segment
  1148.         alpha=((xpr-xp1)*x12+(ypr-yp1)*y12)/sr2
  1149. C
  1150. C       Ecrire la relation si point legerement (1D-4) en dehors du segment
  1151.         pond = 1.d0
  1152.         if (alpha.lt.0.d0) pond = 1.D0 + alpha*1.D4
  1153.         if (alpha.gt.1.d0) pond = 1.D0 + (1.D0 - alpha)*1.D4
  1154.         pond = max(pond,0.D0)
  1155.         pond = min(pond,1.D0)
  1156.         if (pond.le.0.d0) goto 520
  1157. *  verif compatibilite avec la normale au poin impactant
  1158.            if (xms(jp)*xn+yms(jp)*yn.gt. -0.0d0) then
  1159. **        write (6,*) ' impos2 normales incompatibes 1 ',
  1160. **   >  jp,xjeu,dpm,dist
  1161.              goto 520
  1162.            endif
  1163.  1954   continue
  1164. *
  1165. * on a un element ou imposer la relation a ajouter
  1166. *
  1167.         nelri2 = nelri2 + 1
  1168.         nelch = nelch + 1
  1169. *
  1170. *   on ajuste les differents segments
  1171. *
  1172.         if (nelri2.gt.nelrig) then
  1173.           nelrig = nelrig + incnel
  1174.           segadj,xmatri
  1175.           nbelem = nbelem + incnel
  1176.           nbnn = 4
  1177.           segadj,meleme
  1178.           nbnn = 1
  1179.           segadj,ipt7
  1180.           n = n + incnel
  1181.           segadj,mpoval
  1182.         endif
  1183. *
  1184. *   on range dans le meleme
  1185. *
  1186.         num(1,nelri2)=ipt1.num(1,iel1)
  1187.         num(2,nelri2)=ip1
  1188.         num(3,nelri2)=ip2
  1189.         num(4,nelri2)=jp
  1190. *
  1191. *   on remplit le xmatri
  1192. *
  1193. *  lambda
  1194.         re(1,1,nelri2)= 0.d0
  1195. *  ip1
  1196.         re(2,1,nelri2)=  -xn * (1.-alpha) * sr * pond
  1197.         re(3,1,nelri2)=  -yn * (1.-alpha) * sr * pond
  1198. *  ip2
  1199.         re(4,1,nelri2)=  -xn * alpha      * sr * pond
  1200.         re(5,1,nelri2)=  -yn * alpha      * sr * pond
  1201. *  jp
  1202.         re(6,1,nelri2)= xn              * sr   * pond
  1203.         re(7,1,nelri2)= yn              * sr   * pond
  1204. *  on transpose
  1205.         re(1,2,nelri2) = re(2,1,nelri2)
  1206.         re(1,3,nelri2) = re(3,1,nelri2)
  1207.         re(1,4,nelri2) = re(4,1,nelri2)
  1208.         re(1,5,nelri2) = re(5,1,nelri2)
  1209.         re(1,6,nelri2) = re(6,1,nelri2)
  1210.         re(1,7,nelri2) = re(7,1,nelri2)
  1211. *  le reste est nul
  1212. *
  1213. *  remplissage du champoint de depi
  1214. *
  1215.         xjeu1 = x1p*xn + y1p*yn
  1216.         xjeu2 = x2p*xn + y2p*yn
  1217.         xjeu  = (1.-alpha)*xjeu1 + alpha * xjeu2
  1218.         ipt7.num(1,nelch) = ipt1.num(1,iel1)
  1219.         vpocha(nelch,1) = (-xjeu - xjr)*sr * pond
  1220.         vpocha(nelch,2) = sr * pond
  1221.         IF (MELVA2.NE.0) THEN  
  1222.           NEL2 = min (iel1,NELA)
  1223.           VPOCHA(nelch,3) = max(MELVA2.VELCHE(NPTELA,NEL2),0.D0)*sr
  1224.         ENDIF  
  1225. **      write(ioimp,*) ' jeu type 2 ',nelri2,nelch,vpocha(nelch,1),
  1226. **   &                 ipt7.num(1,nelch),ip1,ip2,jp
  1227.  520  continue
  1228.  519  continue
  1229. *
  1230. *      write (ioimp,*) ' nb relation type 2 ',nelri2
  1231. *
  1232. * Ajustement au plus juste puis desactivation des segments lies
  1233. * a la rigidite du type 2
  1234.       if (nelri2.ne.nelrig) then
  1235.         nelrig = nelri2
  1236.         segadj,xmatri
  1237.         nbelem = nelri2
  1238.         nbnn = 4
  1239.         segadj,meleme
  1240.       endif
  1241.       segdes,xmatri
  1242. *
  1243. *  si il n'y a rien on dit que pas unilateral pour pas passer dans unilater
  1244. *     if (nbelem.eq.0) irigel(6,nrigel) = 0
  1245. *
  1246.  700  CONTINUE
  1247. *
  1248.       GOTO 1000
  1249. *
  1250. *----------------------------
  1251. *    on renvoie le resultat
  1252. *----------------------------
  1253.  1000 CONTINUE
  1254.       segsup mfopa2
  1255. *
  1256. * Ajustement au plus juste du chpoint de depi (jeu) : mpoval et ipt7
  1257. * puis desactivation du chpoint
  1258.       if (n.ne.nelch) then
  1259.         n = nelch
  1260.         segadj,mpoval
  1261.         nbnn   = 1
  1262.         nbelem = nelch
  1263.         nbsous = 0
  1264.         nbref  = 0
  1265.         segadj,ipt7
  1266.       endif
  1267. * Desctivation de la matrice de raideur de contact
  1268.       segdes,mrigid
  1269. *
  1270. * Reunion des relations portant sur le meme multiplicateur de lagrange
  1271. *
  1272.       call impofu(MRIGID,MCHPOI)
  1273. C
  1274.  900  continue
  1275.       end
  1276.  
  1277.  
  1278.  

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