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

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