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

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