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

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