Télécharger impo32.eso

Retour à la liste

Numérotation des lignes :

  1. C IMPO32 SOURCE PV 17/03/22 21:15:04 9370
  2.  
  3. * impo bloc en 3D
  4.  
  5. SUBROUTINE IMPO32(ipt6)
  6.  
  7. IMPLICIT INTEGER(I-N)
  8. IMPLICIT REAL*8 (A-H,O-Z)
  9.  
  10. -INC CCOPTIO
  11. -INC CCREEL
  12. -INC CCGEOME
  13.  
  14. -INC SMCOORD
  15. -INC SMELEME
  16. -INC SMRIGID
  17. -INC SMCHPOI
  18. -INC SMCHAML
  19.  
  20. * Contact en formulation forte :
  21. * Segments utiles dans les cas dits pathologiques
  22. segment mfopa1
  23. integer lpoin(nfopa1,3)
  24. real*8 xpoin(nfopa1,3)
  25. real*8 zje1(nfopa1)
  26. endsegment
  27. segment mfopa2
  28. integer lsegm(nfopa2,4)
  29. real*8 xsegm(nfopa2,3)
  30. real*8 zje2(nfopa2)
  31. endsegment
  32. * directions moyennes aux sommets
  33. segment mfopa3
  34. real*8 xms(nbpts),yms(nbpts),zms(nbpts)
  35. endsegment
  36.  
  37. * Contact en formulation faible :
  38. * Segment utile pour le calcul de l'intersection des 2 triangles
  39. * projetes dans le plan de contact intermediaire
  40. segment mfaible
  41. real*8 cPT1C(3,3),cPT2C(3,3), cPT1A(4,2),cPT2A(4,2)
  42. real*8 cPIn0(6,2), cPIn(6,2), test(6)
  43. real*8 vT1A(3,2), vT2A(3,2)
  44. real*8 SuT1A, SuT2A, SuPIn
  45. *DBG-F real*8 xGIn,yGIn, b1T1,b2T1,b3T1, b1T2,b2T2,b3T2
  46. endsegment
  47.  
  48. PARAMETER ( X1s3=0.333333333333333333333333333333333333333333D0 )
  49.  
  50. character*4 modepl(3),moforc(3)
  51.  
  52. data modepl /'UX ','UY ','UZ '/
  53. data moforc /'FX ','FY ','FZ '/
  54. *
  55. idimp1 = IDIM + 1
  56. *
  57. * Lecture du maillage support des conditions de contact
  58. *
  59. call lirobj('MAILLAGE',ipt1,1,iretou)
  60. isous=0
  61. call lirent(isous,0,iretou)
  62. mchel1=0
  63. call lirent(mchel1,0,iretou)
  64. if (ierr.ne.0) return
  65. *
  66. * Activation du maillage et petites verifications
  67. *
  68. segact ipt1
  69. nbno1 = ipt1.num(/1)
  70. nbel1 = ipt1.num(/2)
  71. if (ipt1.lisous(/1).ne.0) call erreur(25)
  72. if (ipt1.itypel.ne.22) call erreur(16)
  73. if (nbno1.ne.5 .and. nbno1.ne.7) call erreur(16)
  74. if (ierr.ne.0) goto 900
  75. *
  76. * Increment sur le nombre d'elements de contact retenus et utilises
  77. * dans le chpoint (mpoval et igeoc) et rigidite pour agrandir les
  78. * longueurs des segments adequats
  79. incnel = 2000
  80. *
  81. * Creation de la raideur des conditions de contact
  82. * Remplissage de l'entete commun
  83. *
  84. nrigel = 1
  85. segini,mrigid
  86. ichole = 0
  87. imgeo1 = 0
  88. imgeo2 = 0
  89. isupeq = 0
  90. iforig = ifomod
  91. coerig(1) = 1.
  92. mtymat='RIGIDITE'
  93. *
  94. * Creation du chpoint de depi
  95. *
  96. nat=1
  97. nsoupo=1
  98. segini mchpoi
  99. mtypoi='DEPIMP'
  100. mochde='engendré par impose'
  101. ifopoi=ifomod
  102. jattri(1)=2
  103. *
  104. nc=1
  105. segini msoupo
  106. ipchp(1)=msoupo
  107. nocomp(1)='FLX '
  108. nocons(1)=' '
  109. *
  110. nbnn =1
  111. nbelem=incnel
  112. nbref =0
  113. nbsous=0
  114. segini ipt7
  115. ipt7.itypel=1
  116. igeoc=ipt7
  117. *
  118. n=incnel
  119. segini mpoval
  120. ipoval=mpoval
  121. *
  122. * Nombre de noeuds dans le chpoint (en totalite) : ipt7 et mpoval
  123. nelch = 0
  124.  
  125. *=======================================================================
  126. * Formulation "forte" (standard) du contact :
  127. *=======================================================================
  128. * Element de contact 3D a 5 noeuds
  129. if (nbno1.eq.7) goto 500
  130. if (mchel1.ne.0) then
  131. segact mchel1
  132. mcham1 = mchel1.ichaml(isous)
  133. segact mcham1
  134. melva1=mcham1.ielval(1)
  135. segact melva1
  136. nptel = melva1.velche(/1)
  137. nel = melva1.velche(/2)
  138. nptel1=min(nptel,4)
  139. endif
  140. * Calcul des directions moyennes
  141. nbpts = xcoor(/1)/(idim+1)
  142. segini mfopa3
  143. ip1o=0
  144. ip2o=0
  145. ip3o=0
  146. do 820 iel=1,nbel1
  147. *
  148. ip1=ipt1.num(2,iel)
  149. ip2=ipt1.num(3,iel)
  150. ip3=ipt1.num(4,iel)
  151. * on suppose que les elements de frottements appuyes sur la même facette se suivent
  152. * et on ne les comptabilise qu'une fois
  153. if (ip1.eq.ip1o.and.ip2.eq.ip2o.and.ip3.eq.ip3o) goto 820
  154. ip1o=ip1
  155. ip2o=ip2
  156. ip3o=ip3
  157. ipv1 = (ip1-1)*idimp1
  158. xp1 = xcoor(ipv1+1)
  159. yp1 = xcoor(ipv1+2)
  160. zp1 = xcoor(ipv1+3)
  161. ipv2 = (ip2-1)*idimp1
  162. xp2 = xcoor(ipv2+1)
  163. yp2 = xcoor(ipv2+2)
  164. zp2 = xcoor(ipv2+3)
  165. ipv3 = (ip3-1)*idimp1
  166. xp3 = xcoor(ipv3+1)
  167. yp3 = xcoor(ipv3+2)
  168. zp3 = xcoor(ipv3+3)
  169. *
  170. * normale au plan (123)
  171. *
  172. x12 = xp2 - xp1
  173. y12 = yp2 - yp1
  174. z12 = zp2 - zp1
  175. x23 = xp3 - xp2
  176. y23 = yp3 - yp2
  177. z23 = zp3 - zp2
  178. xn = y12*z23 - z12*y23
  179. yn = z12*x23 - x12*z23
  180. zn = x12*y23 - y12*x23
  181. sn = xn*xn + yn*yn + zn*zn
  182. sn = max(sqrt(abs(sn)),xpetit*1d10)
  183. xn = xn/sn
  184. yn = yn/sn
  185. zn = zn/sn
  186. xms(ip1)=xms(ip1)+xn
  187. yms(ip1)=yms(ip1)+yn
  188. zms(ip1)=zms(ip1)+zn
  189. xms(ip2)=xms(ip2)+xn
  190. yms(ip2)=yms(ip2)+yn
  191. zms(ip2)=zms(ip2)+zn
  192. xms(ip3)=xms(ip3)+xn
  193. yms(ip3)=yms(ip3)+yn
  194. zms(ip3)=zms(ip3)+zn
  195. 820 continue
  196. do 821 ip = 1,nbpts
  197. sn = xms(ip)*xms(ip)+yms(ip)*yms(ip)+zms(ip)*zms(ip)
  198. sn = max(sqrt(abs(sn)),xpetit*1d10)
  199. xms(ip)=xms(ip)/sn
  200. yms(ip)=yms(ip)/sn
  201. zms(ip)=zms(ip)/sn
  202. 821 continue
  203. *
  204. * Nombre d'iterations pour la detection de la direction du contact
  205. nbiter=6
  206. *
  207. * Relation du type 3 : noeud-triangle
  208. *-------------------------------------
  209. *
  210. * creation du meleme associe a la relation
  211. *
  212. nbsous=0
  213. nbref =0
  214. nbnn =5
  215. nbelem=incnel
  216. segini meleme
  217. itypel=22
  218. irigel(1,nrigel)=meleme
  219. *
  220. * creation du descriptif commun a toutes les raideurs
  221. *
  222. nligrp = 13
  223. nligrd = nligrp
  224. segini,descr
  225. lisinc(1)='LX '
  226. lisdua(1)='FLX '
  227. noelep(1)=1
  228. noeled(1)=1
  229. do 10 i=2, nligrp, 3
  230. lisinc(i) =modepl(1)
  231. lisinc(i+1)=modepl(2)
  232. lisinc(i+2)=modepl(3)
  233. lisdua(i) =moforc(1)
  234. lisdua(i+1)=moforc(2)
  235. lisdua(i+2)=moforc(3)
  236. noelep(i) =(i+4)/3
  237. noelep(i+1)=noelep(i)
  238. noelep(i+2)=noelep(i)
  239. noeled(i) =noelep(i)
  240. noeled(i+1)=noelep(i)
  241. noeled(i+2)=noelep(i)
  242. 10 continue
  243. segdes,descr
  244. irigel(3,nrigel) = descr
  245. *
  246. * creation du segment xmatri
  247. *
  248. nelrig = incnel
  249. segini,xmatri
  250. irigel(4,nrigel) = xmatri
  251. *
  252. * ce qu'on cree est unilateral
  253. *
  254. irigel(6,nrigel)=1
  255. *
  256. * ce qu'on cree est symetrique
  257. *
  258. irigel(7,nrigel)=0
  259. *
  260. * Nombre d'elements dans la rigidite du type 3
  261. nelri3 = 0
  262.  
  263. * Boucle sur le maillage support du contact(frottement)
  264. *
  265. DO 20 iel=1,nbel1
  266. *
  267. xjr = 0d0
  268. ip1=ipt1.num(2,iel)
  269. ip2=ipt1.num(3,iel)
  270. ip3=ipt1.num(4,iel)
  271. * Recuperation des coordonees des noeuds du triangle
  272. ipv = (ip1-1)*idimp1
  273. xp1 = xcoor(ipv+1)
  274. yp1 = xcoor(ipv+2)
  275. zp1 = xcoor(ipv+3)
  276. ipv = (ip2-1)*idimp1
  277. xp2 = xcoor(ipv+1)
  278. yp2 = xcoor(ipv+2)
  279. zp2 = xcoor(ipv+3)
  280. ipv = (ip3-1)*idimp1
  281. xp3 = xcoor(ipv+1)
  282. yp3 = xcoor(ipv+2)
  283. zp3 = xcoor(ipv+3)
  284. * Centre de gravite du triangle
  285. xg = (xp1+xp2+xp3) /3D0
  286. yg = (yp1+yp2+yp3) /3d0
  287. zg = (zp1+zp2+zp3) /3d0
  288. jp =ipt1.num(5,iel)
  289. ipv = (jp-1)*idimp1
  290. xp = xcoor(ipv+1)
  291. yp = xcoor(ipv+2)
  292. zp = xcoor(ipv+3)
  293. if (mchel1.ne.0) then
  294. nel1 = min (iel,nel)
  295. xjr = melva1.velche(nptel1,nel1)
  296. endif
  297. * Verification que pas relation du point sur lui meme
  298. if (jp.eq.ip1) goto 20
  299. if (jp.eq.ip2) goto 20
  300. if (jp.eq.ip3) goto 20
  301. *
  302. * Verification si autre point dans la zone de selection
  303. * verif distance au centre de gravite
  304. d1 = (xg-xp1)**2 + (yg-yp1)**2 + (zg-zp1)**2
  305. d2 = (xg-xp2)**2 + (yg-yp2)**2 + (zg-zp2)**2
  306. d3 = (xg-xp3)**2 + (yg-yp3)**2 + (zg-zp3)**2
  307. dist = max(d1,d2,d3)
  308. dp = ((xg-xp)**2 + (yg-yp)**2 + (zg-zp)**2)
  309. if (dp.gt.dist*16) goto 20
  310. C*DBG write(ioimp,*) 'contact test distance ok',dp,d1,d2,d3
  311. dist = SQRT(dist)
  312. * verif position par rapport aux cotes
  313. * cote 1 2
  314. x12 = xp2 - xp1
  315. y12 = yp2 - yp1
  316. z12 = zp2 - zp1
  317. dv = sqrt((x12*x12)+(y12*y12)+(z12*z12))
  318. if (.not.(dv.lt.1d-10) .and. .not.(dv.gt.1d-10))
  319. & write(ioimp,*) 'Prob 1.1 - impo32'
  320. if (abs(dv).le.xpetit) write(ioimp,*) 'Prob 1.2 - impo32'
  321. xv = x12 / dv
  322. yv = y12 / dv
  323. zv = z12 / dv
  324. scal = (xp-xp1)*xv + (yp-yp1)*yv + (zp-zp1)*zv
  325. xm = xp1 + scal*xv
  326. ym = yp1 + scal*yv
  327. zm = zp1 + scal*zv
  328. scal = (xg-xp1)*xv + (yg-yp1)*yv + (zg-zp1)*zv
  329. xn = (xg-xp1) - scal*xv
  330. yn = (yg-yp1) - scal*yv
  331. zn = (zg-zp1) - scal*zv
  332. dn = sqrt(xn*xn+yn*yn+zn*zn)
  333. dpm = sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2)
  334. scal = (xp-xm)*xn + (yp-ym)*yn + (zp-zm)*zn
  335. if (dpm.gt.1d-9*dv .and. scal.lt.-0.5d0*dn*dpm) goto 20
  336. dpm=max(xpetit,dpm)
  337. * write(ioimp,*) 'contact test position 1 ok',
  338. * & scal/(dn*dpm),scal,dn,dpm
  339. * cote 2 3
  340. x23 = xp3 - xp2
  341. y23 = yp3 - yp2
  342. z23 = zp3 - zp2
  343. dv = sqrt((x23*x23)+(y23*y23)+(z23*z23))
  344. if (.not.(dv.lt.1d-10) .and. .not.(dv.gt.1d-10))
  345. & write(ioimp,*) 'Prob 2.1 - impo32'
  346. if (abs(dv).le.xpetit) write(ioimp,*) 'Prob 2.2 - impo32'
  347. xv = x23 / dv
  348. yv = y23 / dv
  349. zv = z23 / dv
  350. scal = (xp-xp2)*xv + (yp-yp2)*yv + (zp-zp2)*zv
  351. xm = xp2 + scal*xv
  352. ym = yp2 + scal*yv
  353. zm = zp2 + scal*zv
  354. scal = (xg-xp2)*xv + (yg-yp2)*yv + (zg-zp2)*zv
  355. xn = (xg-xp2) - scal*xv
  356. yn = (yg-yp2) - scal*yv
  357. zn = (zg-zp2) - scal*zv
  358. dn = sqrt(xn*xn+yn*yn+zn*zn)
  359. dpm = sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2)
  360. scal = (xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn
  361. if (dpm.gt.1d-9*dv .and. scal.lt.-0.5d0*dn*dpm) goto 20
  362. dpm=max(xpetit,dpm)
  363. * write(ioimp,*) 'contact test position 2 ok',
  364. * & scal/(dn*dpm),scal,dn,dpm
  365. * cote 3 1
  366. x31 = xp1 - xp3
  367. y31 = yp1 - yp3
  368. z31 = zp1 - zp3
  369. dv = sqrt((x31*x31)+(y31*y31)+(z31*z31))
  370. if (.not.(dv.lt.1d-10) .and. .not.(dv.gt.1d-10))
  371. & write(ioimp,*) 'Prob 3.1 - impo32'
  372. if (abs(dv).le.xpetit) write(ioimp,*) 'Prob 3.2 - impo32'
  373. xv = x31 / dv
  374. yv = y31 / dv
  375. zv = z31 / dv
  376. scal = (xp-xp3)*xv + (yp-yp3)*yv + (zp-zp3)*zv
  377. xm = xp3 + scal*xv
  378. ym = yp3 + scal*yv
  379. zm = zp3 + scal*zv
  380. scal = (xg-xp3)*xv + (yg-yp3)*yv + (zg-zp3)*zv
  381. xn = (xg-xp3) - scal*xv
  382. yn = (yg-yp3) - scal*yv
  383. zn = (zg-zp3) - scal*zv
  384. dn = sqrt(xn*xn+yn*yn+zn*zn)
  385. dpm = sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2)
  386. scal = (xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn
  387. if (dpm.gt.1d-9*dv .and .scal.lt.-0.5d0*dn*dpm) goto 20
  388. dpm=max(xpetit,dpm)
  389. * write(ioimp,*) 'contact test position 3 ok',
  390. * & scal/(dn*dpm),scal,dn,dpm
  391. *
  392. * on a un point ou imposer la relation
  393. C*DBG write(ioimp,*) 'contact potentiel '
  394. * Quelle est la relation ???
  395. *
  396. * direction de la relation = normale au plan (123)
  397. *
  398. xn = y12*z23 - z12*y23
  399. yn = z12*x23 - x12*z23
  400. zn = x12*y23 - y12*x23
  401. dn = sqrt(xn*xn+yn*yn+zn*zn)
  402. if (.not.(dn.lt.1d-10).and. .not.(dn.gt.1d-10))
  403. & write(ioimp,*) 'Prob 4.1 - impo32'
  404. if (abs(dn).le.xpetit) write(ioimp,*) 'Prob 4.2 - impo32'
  405. xn = xn / dn
  406. yn = yn / dn
  407. zn = zn / dn
  408. *
  409. * Distance (jeu) du point jp au plan de contact
  410. * Puis calcul de la projection sur le plan de contact
  411. *
  412. xjeu = (xp-xg)*xn + (yp-yg)*yn + (zp-zg)*zn - xjr
  413. * verif bon cote (a la tolerance pres)
  414. * write (6,*) ' xjeu ',xjeu,dist,iel
  415.  
  416. if (xjeu.lt.-3*dist) goto 20
  417. if (xjeu.gt.3*dist) goto 20
  418.  
  419. *
  420. * Recherche de ses coordonnées barycentriques
  421. * qui sont les surfaces des sous triangles
  422. *
  423. xb1 = (yp-yp2)*(zp-zp3)-(zp-zp2)*(yp-yp3)
  424. yb1 = (zp-zp2)*(xp-xp3)-(xp-xp2)*(zp-zp3)
  425. zb1 = (xp-xp2)*(yp-yp3)-(yp-yp2)*(xp-xp3)
  426. xb2 = (yp-yp3)*(zp-zp1)-(zp-zp3)*(yp-yp1)
  427. yb2 = (zp-zp3)*(xp-xp1)-(xp-xp3)*(zp-zp1)
  428. zb2 = (xp-xp3)*(yp-yp1)-(yp-yp3)*(xp-xp1)
  429. xb3 = (yp-yp1)*(zp-zp2)-(zp-zp1)*(yp-yp2)
  430. yb3 = (zp-zp1)*(xp-xp2)-(xp-xp1)*(zp-zp2)
  431. zb3 = (xp-xp1)*(yp-yp2)-(yp-yp1)*(xp-xp2)
  432. iter = 0
  433. 21 continue
  434. b1 = xb1*xn + yb1*yn + zb1*zn
  435. b2 = xb2*xn + yb2*yn + zb2*zn
  436. b3 = xb3*xn + yb3*yn + zb3*zn
  437. bt = b1 + b2 + b3
  438. bsurf = bt / 2
  439. if (.not.(bt.lt.1d-10) .and. .not.(bt.gt.1d-10))
  440. & write(ioimp,*) 'Prob 5.1 - impo32'
  441. if (abs(bt).le.xpetit) write(ioimp,*) 'Prob 5.2 - impo32'
  442. b1 = b1 / bt
  443. b2 = b2 / bt
  444. b3 = b3 / bt
  445. * recalcul de la direction et du jeu en fonction des normales aux sommets
  446. xn = b1*xms(ip1)+b2*xms(ip2)+b3*xms(ip3)
  447. yn = b1*yms(ip1)+b2*yms(ip2)+b3*yms(ip3)
  448. zn = b1*zms(ip1)+b2*zms(ip2)+b3*zms(ip3)
  449. sn = sqrt (xn*xn + yn*yn + zn*zn)
  450. xn=xn/sn
  451. yn=yn/sn
  452. zn=zn/sn
  453. * recalcul en fonction de la nouvelle normale
  454. iter=iter+1
  455. if (iter.le.nbiter) goto 21
  456. * verif compatibilite avec la normale au poin impactant
  457. if (xms(jp)*xn+yms(jp)*yn+zms(jp)*zn.gt. 0.0d0) then
  458. ** write (6,*) ' impo32 normales incompatibes 1 '
  459. goto 20
  460. endif
  461. xjeu1 = (xp-xp1)*xn + (yp-yp1)*yn + (zp-zp1)*zn
  462. xjeu2 = (xp-xp2)*xn + (yp-yp2)*yn + (zp-zp2)*zn
  463. xjeu3 = (xp-xp3)*xn + (yp-yp3)*yn + (zp-zp3)*zn
  464. xjeu = xjeu1 * b1 + xjeu2 * b2 + xjeu3 * b3 - xjr
  465. ** write (6,*) ' normale ',xn,yn,zn,' jeu ',xjeu1,xjeu2,xjeu3
  466.  
  467. C*DBG write(ioimp,*) ' b1 b2 b3 ',b1,b2,b3
  468. * creation du support du multiplicateur de lagrange si necessaire
  469. nbpts = xcoor(/1)/(idim+1)
  470. if (ipt1.num(1,iel).eq.ilgni) then
  471. nbpts = nbpts + 1
  472. segadj mcoord
  473. segact ipt1*mod
  474. ipt1.num(1,iel)=nbpts
  475. if (ipt6.ne.0) ipt6.num(1,iel)=nbpts
  476. endif
  477. segact mcoord*mod
  478. xcoor(nbpts*(idim+1)) = dist
  479. 1954 continue
  480. * test qu'on n'est pas trop a l'exterieur
  481. bm=1.d0
  482. if (b1.lt.0.d0) bm=bm+b1*1.0d0
  483. if (b2.lt.0.d0) bm=bm+b2*1.0d0
  484. if (b3.lt.0.d0) bm=bm+b3*1.0d0
  485. bsurf = bsurf * bm
  486. bb1=max(b1,0.d0)
  487. bb2=max(b2,0.d0)
  488. bb3=max(b3,0.d0)
  489.  
  490. if (xjeu.gt.0.d0)then
  491. xjeu = xjeu1 * abs(bb1) + xjeu2 * abs(bb2) + xjeu3 * abs(bb3)
  492. > -xjr
  493. endif
  494. if (bsurf.le.0.d0) goto 20
  495. *
  496. * Ajout d'une relation noeud-triangle
  497. *
  498. nelri3 = nelri3 + 1
  499. nelch = nelch + 1
  500. *
  501. * on ajuste les differents segments si necesaire
  502. *
  503. if (nelri3.gt.nelrig) then
  504. nelrig = nelrig + incnel
  505. segadj,xmatri
  506. nbelem = nbelem + incnel
  507. nbnn = 5
  508. segadj,meleme
  509. nbnn = 1
  510. segadj,ipt7
  511. n = n + incnel
  512. segadj,mpoval
  513. endif
  514. *
  515. * Mise a jour du meleme
  516. num(1,nelri3) = ipt1.num(1,iel)
  517. num(2,nelri3) = ip1
  518. num(3,nelri3) = ip2
  519. num(4,nelri3) = ip3
  520. num(5,nelri3) = jp
  521. *
  522. * Mise a jour de xmatri
  523. * lambda
  524. re( 1,1,nelri3) = 0d0
  525. * ip1
  526. re( 2,1,nelri3) = xn * b1 * bsurf
  527. re( 3,1,nelri3) = yn * b1 * bsurf
  528. re( 4,1,nelri3) = zn * b1 * bsurf
  529. * ip2
  530. re( 5,1,nelri3) = xn * b2 * bsurf
  531. re( 6,1,nelri3) = yn * b2 * bsurf
  532. re( 7,1,nelri3) = zn * b2 * bsurf
  533. * ip3
  534. re( 8,1,nelri3) = xn * b3 * bsurf
  535. re( 9,1,nelri3) = yn * b3 * bsurf
  536. re(10,1,nelri3) = zn * b3 * bsurf
  537. * jp
  538. re(11,1,nelri3) = -xn * bsurf
  539. re(12,1,nelri3) = -yn * bsurf
  540. re(13,1,nelri3) = -zn * bsurf
  541. * on transpose
  542. do 30 ic = 2, nligrp
  543. re(1,ic,nelri3) = re(ic,1,nelri3)
  544. 30 continue
  545. * le reste est nul
  546. *
  547. * remplissage du champoint de depi (jeu)
  548. *
  549. ipt7.num(1,nelch) = ipt1.num(1,iel)
  550. vpocha(nelch,1) = xjeu * bsurf
  551. * write(ioimp,*) ' jeu ',xjeu,ip1,ip2,ip3,jp,b1,b2,b3,xn,yn,zn
  552. *
  553. 20 CONTINUE
  554. *
  555. * Ajustement au plus juste puis desactivation des segments lies
  556. * a la rigidite du type 3
  557. if (nelri3.ne.nelrig) then
  558. nelrig = nelri3
  559. segadj,xmatri
  560. nbelem = nelri3
  561. nbnn = 5
  562. segadj,meleme
  563. endif
  564. segdes,meleme,xmatri
  565. *
  566. c write(ioimp,*) ' nb relation type 3 ',nelri3
  567. C*DBG if (nelri3.eq.0) irigel(6,nrigel)=0
  568. *
  569. * Fin du traitement de la formulation standard/forte
  570. * ----------------------------------------------------
  571. 300 CONTINUE
  572. * Destruction des segments locaux devenus inutiles
  573. segsup mfopa3
  574.  
  575. GOTO 1000
  576.  
  577. *=======================================================================
  578. *= Formulation "faible" 3D : relation triangle-triangle
  579. *=======================================================================
  580. * Element de contact 3D a 7 noeuds (1+6)
  581. 500 CONTINUE
  582. *
  583. * Petit segment de travail
  584. segini,mfaible
  585. *
  586. * Relation du type 0 : triangle-triangle (faible)
  587. *---------------------------------------
  588. *
  589. * Creation du meleme associe a la relation
  590. *
  591. nbnn = 7
  592. nbelem = incnel
  593. nbsous = 0
  594. nbref = 0
  595. segini,meleme
  596. itypel = 22
  597. irigel(1,nrigel) = meleme
  598. *
  599. * Creation du descriptif commun a toutes les raideurs
  600. *
  601. nligrp = 19
  602. nligrd = nligrp
  603. segini,descr
  604. lisinc(1) = 'LX '
  605. lisdua(1) = 'FLX '
  606. noelep(1) = 1
  607. noeled(1) = 1
  608. do 510 i = 2, nligrp, 3
  609. lisinc(i ) = modepl(1)
  610. lisinc(i+1) = modepl(2)
  611. lisinc(i+2) = modepl(3)
  612. lisdua(i ) = moforc(1)
  613. lisdua(i+1) = moforc(2)
  614. lisdua(i+2) = moforc(3)
  615. noelep(i ) = (i+5)/3
  616. noelep(i+1) = noelep(i)
  617. noelep(i+2) = noelep(i)
  618. noeled(i ) = noelep(i)
  619. noeled(i+1) = noelep(i)
  620. noeled(i+2) = noelep(i)
  621. 510 continue
  622. segdes,descr
  623. irigel(3,nrigel) = descr
  624. *
  625. * creation du segment xmatri
  626. *
  627. nelrig = incnel
  628. segini,xmatri
  629. irigel(4,nrigel) = xmatri
  630. *
  631. * ce qu'on cree est unilateral
  632. *
  633. irigel(6,nrigel) = 1
  634. *
  635. * ce qu'on cree est symetrique
  636. *
  637. irigel(7,nrigel) = 0
  638. *
  639. * Nombre d'elements crees dans meleme=irigel(1,nrigel), ipt7 et mpoval
  640. nelri0 = 0
  641. if (mchel1.ne.0) then
  642. segact mchel1
  643. mcham1 = mchel1.ichaml(isous)
  644. segact mcham1
  645. melva1=mcham1.ielval(1)
  646. segact melva1
  647. nptel = melva1.velche(/1)
  648. nel = melva1.velche(/2)
  649. nptel1=min(nptel,4)
  650. endif
  651. *
  652. * Boucle sur les elements du maillage de contact/frottement "faible"
  653. *
  654. DO 520 iel = 1, nbel1
  655. *
  656. xjr = 0d0
  657. ip1 = ipt1.num(2,iel)
  658. ip2 = ipt1.num(3,iel)
  659. ip3 = ipt1.num(4,iel)
  660. jp1 = ipt1.num(5,iel)
  661. jp2 = ipt1.num(6,iel)
  662. jp3 = ipt1.num(7,iel)
  663. if (mchel1.ne.0) then
  664. nel1 = min (iel,nel)
  665. xjr = melva1.velche(nptel1,nel1)
  666. endif
  667. C*DBG write(ioimp,*) iel,ip1,ip2,ip3,jp1,jp2,jp3,ipt1.num(1,iel)
  668. *
  669. * Verification (provisoire) qu'il n'y a pas de noeud double
  670. if (ip1.eq.jp1) goto 520
  671. if (ip1.eq.jp2) goto 520
  672. if (ip1.eq.jp3) goto 520
  673. if (ip2.eq.jp1) goto 520
  674. if (ip2.eq.jp2) goto 520
  675. if (ip2.eq.jp3) goto 520
  676. if (ip3.eq.jp1) goto 520
  677. if (ip3.eq.jp2) goto 520
  678. if (ip3.eq.jp3) goto 520
  679. *
  680. * Definition du triangle T1(ip1,ip2,ip3)
  681. * - Coordonnees des noeuds de T1
  682. ipv = (ip1-1)*idimp1
  683. xp1 = xcoor(ipv+1)
  684. yp1 = xcoor(ipv+2)
  685. zp1 = xcoor(ipv+3)
  686. ipv = (ip2-1)*idimp1
  687. xp2 = xcoor(ipv+1)
  688. yp2 = xcoor(ipv+2)
  689. zp2 = xcoor(ipv+3)
  690. ipv = (ip3-1)*idimp1
  691. xp3 = xcoor(ipv+1)
  692. yp3 = xcoor(ipv+2)
  693. zp3 = xcoor(ipv+3)
  694. * - Barycentre de T1
  695. xgT1 = (xp1 + xp2 + xp3) * X1s3
  696. ygT1 = (yp1 + yp2 + yp3) * X1s3
  697. zgT1 = (zp1 + zp2 + zp3) * X1s3
  698. * - Normale au triangle T1
  699. xnT1 = (yp1-yp2)*(zp2-zp3) - (zp1-zp2)*(yp2-yp3)
  700. ynT1 = (zp1-zp2)*(xp2-xp3) - (xp1-xp2)*(zp2-zp3)
  701. znT1 = (xp1-xp2)*(yp2-yp3) - (yp1-yp2)*(xp2-xp3)
  702. dnT1 = SQRT(xnT1*xnT1+ynT1*ynT1+znT1*znT1)
  703. IF (.not.(dnT1.lt.1d-10) .and. .not.(dnT1.gt.1d-10))
  704. & write(ioimp,*) 'FAIBle - 1.1 - impo32'
  705. IF (abs(dnT1).le.xpetit) write(ioimp,*) 'FAIBle - 1.2 - impo32'
  706. xnT1 = xnT1 / dnT1
  707. ynT1 = ynT1 / dnT1
  708. znT1 = znT1 / dnT1
  709. C*DBG write(ioimp,*) 'Normale au plan T1',xnT1,ynT1,znT1,dnT1
  710. * - Rayon de la sphere "englobante" de l'element (centre=barycentre)
  711. d1 = (xgT1-xp1)**2 + (ygT1-yp1)**2 + (zgT1-zp1)**2
  712. d2 = (xgT1-xp2)**2 + (ygT1-yp2)**2 + (zgT1-zp2)**2
  713. d3 = (xgT1-xp3)**2 + (ygT1-yp3)**2 + (zgT1-zp3)**2
  714. RayT1 = SQRT(max(d1,d2,d3))
  715. *
  716. * Definition du triangle T2(jp1,jp2,jp3)
  717. * - Coordonnees des noeuds de T2
  718. ipv = (jp1-1)*idimp1
  719. xq1 = xcoor(ipv+1)
  720. yq1 = xcoor(ipv+2)
  721. zq1 = xcoor(ipv+3)
  722. ipv = (jp2-1)*idimp1
  723. xq2 = xcoor(ipv+1)
  724. yq2 = xcoor(ipv+2)
  725. zq2 = xcoor(ipv+3)
  726. ipv = (jp3-1)*idimp1
  727. xq3 = xcoor(ipv+1)
  728. yq3 = xcoor(ipv+2)
  729. zq3 = xcoor(ipv+3)
  730. * - Barycentre de T2
  731. xgT2 = (xq1 + xq2 + xq3) * X1s3
  732. ygT2 = (yq1 + yq2 + yq3) * X1s3
  733. zgT2 = (zq1 + zq2 + zq3) * X1s3
  734. * - Normale au triangle T2
  735. xnT2 = (yq1-yq2)*(zq2-zq3) - (zq1-zq2)*(yq2-yq3)
  736. ynT2 = (zq1-zq2)*(xq2-xq3) - (xq1-xq2)*(zq2-zq3)
  737. znT2 = (xq1-xq2)*(yq2-yq3) - (yq1-yq2)*(xq2-xq3)
  738. dnT2 = SQRT(xnT2*xnT2+ynT2*ynT2+znT2*znT2)
  739. IF (.not.(dnT2.lt.1d-10) .and. .not.(dnT2.gt.1d-10))
  740. & write(ioimp,*) 'FAIBle - 2.1 - impo32'
  741. IF (abs(dnT2).le.xpetit) write(ioimp,*) 'FAIBle - 2.2 - impo32'
  742. xnT2 = xnT2 / dnT2
  743. ynT2 = ynT2 / dnT2
  744. znT2 = znT2 / dnT2
  745. C*DBG write(ioimp,*) 'Normale au plan T2',xnT2,ynT2,znT2,dnT2
  746. * - Rayon de la sphere "englobante" de l'element (centre=barycentre)
  747. d1 = (xgT2-xq1)**2 + (ygT2-yq1)**2 + (zgT2-zq1)**2
  748. d2 = (xgT2-xq2)**2 + (ygT2-yq2)**2 + (zgT2-zq2)**2
  749. d3 = (xgT2-xq3)**2 + (ygT2-yq3)**2 + (zgT2-zq3)**2
  750. RayT2 = SQRT(max(d1,d2,d3))
  751.  
  752. * Orientation respective correcte des 2 normales
  753. * Les triangles T1 et T2 etant decrits dans le meme sens ("prisme"),
  754. * le produit scalaire de leur normale doit etre positif ou nul !
  755. scal = xnT1 * xnT2 + ynT1 * ynT2 + znT1 * znT2
  756. C*DBG write(ioimp,*) 'Prod.scal des normales',scal
  757. if (scal.lt.0.) goto 520
  758. *
  759. * Distance entre les centres de gravite
  760. * Les spheres "englobantes" doivent s'intersecter pour contact potentiel
  761. dist = SQRT((xgT2-xgT1)**2 + (ygT2-ygT1)**2 + (zgT2-zgT1)**2)
  762. IF (dist*4.GT.RayT1+RayT2) goto 520
  763. *
  764. * Critere d'acceptation de l'element de contact
  765. *
  766. * Definition du plan de contact (surface intermediaire)
  767. * - Normale (unitaire) du plan de contact
  768. xnC = xnT1 + xnT2
  769. ynC = ynT1 + ynT2
  770. znC = znT1 + znT2
  771. dnC = SQRT(xnC*xnC + ynC*ynC + znC*znC)
  772. IF (.not.(dnC.lt.1d-10) .and. .not.(dnC.GT.1D-10))
  773. & write(ioimp,*) 'FAIBle - 3.1 - impo32'
  774. IF (abs(dnC).le.xpetit) write(ioimp,*) 'FAIBle - 3.2 - impo32'
  775. xnC = xnC / dnC
  776. ynC = ynC / dnC
  777. znC = znC / dnC
  778. C*DBG write(ioimp,*) 'Normale au plan Contact',xnC,ynC,znC,dnC
  779.  
  780. * - Point definissant le plan de contact
  781. * = milieu des barycentres des triangles T1 et T2
  782. xgC = (xgT1 + xgT2) * 0.5
  783. ygC = (ygT1 + ygT2) * 0.5
  784. zgC = (zgT1 + zgT2) * 0.5
  785. * - "Distance" de reference au plan de contact
  786. dgC = xgC*xnC + ygC*ynC + zgC*znC
  787. *
  788. * Distance "signee" des triangles T1 et T2 au plan de contact
  789. dp1C = (xp1*xnC + yp1*ynC + zp1*znC) - dgC
  790. dp2C = (xp2*xnC + yp2*ynC + zp2*znC) - dgC
  791. dp3C = (xp3*xnC + yp3*ynC + zp3*znC) - dgC
  792. dq1C = (xq1*xnC + yq1*ynC + zq1*znC) - dgC
  793. dq2C = (xq2*xnC + yq2*ynC + zq2*znC) - dgC
  794. dq3C = (xq3*xnC + yq3*ynC + zq3*znC) - dgC
  795. * Projection des triangles T1 et T2 sur le plan de contact
  796. * T1C = triangle T1 projete sur plan de Contact
  797. cPT1C(1,1) = xp1 - dp1C * xnC
  798. cPT1C(1,2) = yp1 - dp1C * ynC
  799. cPT1C(1,3) = zp1 - dp1C * znC
  800. cPT1C(2,1) = xp2 - dp2C * xnC
  801. cPT1C(2,2) = yp2 - dp2C * ynC
  802. cPT1C(2,3) = zp2 - dp2C * znC
  803. cPT1C(3,1) = xp3 - dp3C * xnC
  804. cPT1C(3,2) = yp3 - dp3C * ynC
  805. cPT1C(3,3) = zp3 - dp3C * znC
  806. * T2C = triangle T2 projete sur plan de Contact
  807. cPT2C(1,1) = xq1 - dq1C * xnC
  808. cPT2C(1,2) = yq1 - dq1C * ynC
  809. cPT2C(1,3) = zq1 - dq1C * znC
  810. cPT2C(2,1) = xq2 - dq2C * xnC
  811. cPT2C(2,2) = yq2 - dq2C * ynC
  812. cPT2C(2,3) = zq2 - dq2C * znC
  813. cPT2C(3,1) = xq3 - dq3C * xnC
  814. cPT2C(3,2) = yq3 - dq3C * ynC
  815. cPT2C(3,3) = zq3 - dq3C * znC
  816. *
  817. * On determine quelle est la composante maximale de la normale pour
  818. * projeter les triangles T1C et T2C sur un plan parallele aux axes
  819. * ("le plus orthogonal" a cette normale), ce qui maximise leur aire.
  820. r_x = abs(xnC)
  821. r_y = abs(ynC)
  822. r_z = abs(znC)
  823. if (r_x.gt.r_y) then
  824. inC = 1
  825. if (r_x.lt.r_z) inC = 3
  826. else
  827. inC = 2
  828. if (r_y.lt.r_z) inC = 3
  829. endif
  830. * Projection des triangles T1C et T2C sur ce plan // aux axes (A)
  831. * T1A = triangle T1C projete sur ce plan (A)
  832. * T2A = triangle T2C projete sur ce plan (A)
  833. * On prend soin d'orienter les triangles T1A et T2A dans le sens
  834. * direct (en tenant compte du signe de la composante normale !)
  835. GOTO (531,532,533), inC
  836. call erreur(5)
  837. return
  838. * Normale NC selon "x" : A =(y,z)
  839. 531 if (xnC.lt.0.) inC = -inC
  840. inX = 2
  841. inY = 3
  842. goto 534
  843. * Normale NC selon "y" : A =(x,z)
  844. 532 if (ynC.lt.0.) inC = -inC
  845. inX = 1
  846. inY = 3
  847. goto 534
  848. * Normale NC selon "z" : A =(x,y)
  849. 533 if (znC.lt.0.) inC = -inC
  850. inX = 1
  851. inY = 2
  852. goto 534
  853. 534 continue
  854. cPT1A(1,1) = cPT1C(1,inX)
  855. cPT1A(1,2) = cPT1C(1,inY)
  856. cPT2A(1,1) = cPT2C(1,inX)
  857. cPT2A(1,2) = cPT2C(1,inY)
  858. if (inC.gt.0) then
  859. cPT1A(2,1) = cPT1C(2,inX)
  860. cPT1A(2,2) = cPT1C(2,inY)
  861. cPT2A(2,1) = cPT2C(2,inX)
  862. cPT2A(2,2) = cPT2C(2,inY)
  863. cPT1A(3,1) = cPT1C(3,inX)
  864. cPT1A(3,2) = cPT1C(3,inY)
  865. cPT2A(3,1) = cPT2C(3,inX)
  866. cPT2A(3,2) = cPT2C(3,inY)
  867. else
  868. cPT1A(2,1) = cPT1C(3,inX)
  869. cPT1A(2,2) = cPT1C(3,inY)
  870. cPT2A(2,1) = cPT2C(3,inX)
  871. cPT2A(2,2) = cPT2C(3,inY)
  872. cPT1A(3,1) = cPT1C(2,inX)
  873. cPT1A(3,2) = cPT1C(2,inY)
  874. cPT2A(3,1) = cPT2C(2,inX)
  875. cPT2A(3,2) = cPT2C(2,inY)
  876. endif
  877. cPT1A(4,1) = cPT1A(1,1)
  878. cPT1A(4,2) = cPT1A(1,2)
  879. * Calcul des cotes des triangles T1A et T2A (12,23,31)
  880. * Utile pour connaitre ensuite la normale aux cotes des triangles
  881. vT1A(1,1) = cPT1A(2,1) - cPT1A(1,1)
  882. vT1A(1,2) = cPT1A(2,2) - cPT1A(1,2)
  883. vT1A(2,1) = cPT1A(3,1) - cPT1A(2,1)
  884. vT1A(2,2) = cPT1A(3,2) - cPT1A(2,2)
  885. vT1A(3,1) = cPT1A(1,1) - cPT1A(3,1)
  886. vT1A(3,2) = cPT1A(1,2) - cPT1A(3,2)
  887. C*nu vT2A(1,1) = cPT2A(2,1) - cPT2A(1,1)
  888. vT2A(1,2) = cPT2A(2,2) - cPT2A(1,2)
  889. C*nu vT2A(2,1) = cPT2A(3,1) - cPT2A(2,1)
  890. vT2A(2,2) = cPT2A(3,2) - cPT2A(2,2)
  891. C*nu vT2A(3,1) = cPT2A(1,1) - cPT2A(3,1)
  892. vT2A(3,2) = cPT2A(1,2) - cPT2A(3,2)
  893. * Calcul de la surface de chacun des triangles T1A et T2A
  894. * (en fait, on calcule le double, mais par la suite on ne considere que
  895. * des rapports de surfaces de triangle)
  896. SuT1A = (cPT1A(2,1)+cPT1A(1,1)) * vT1A(1,2)
  897. & + (cPT1A(3,1)+cPT1A(2,1)) * vT1A(2,2)
  898. & + (cPT1A(1,1)+cPT1A(3,1)) * vT1A(3,2)
  899. SuT2A = (cPT2A(2,1)+cPT2A(1,1)) * vT2A(1,2)
  900. & + (cPT2A(3,1)+cPT2A(2,1)) * vT2A(2,2)
  901. & + (cPT2A(1,1)+cPT2A(3,1)) * vT2A(3,2)
  902. C*DBG write(ioimp,*) 'Surfaces T1A - T2A :', 0.5*SuT1A,0.5*SuT2A
  903. C*DBG if (SuT1A.gt.100.*SuT2A .or. SuT2A.gt.100.*SuT1A)
  904. C*DBG & write(ioimp,*) 'Rapport des surfaces tres important !'
  905.  
  906. * On initialise l'intersection des triangles avec T2A
  907. nPIn0 = 3
  908. do 540 i = 1, nPIn0
  909. cPIn0(i,1) = cPT2A(i,1)
  910. cPIn0(i,2) = cPT2A(i,2)
  911. 540 continue
  912. * Determination de l'intersection des triangles T1A et T2A
  913. * en regardant progressivement l'intersection avec les cotes de T1A
  914. * Note : On utilise le fait que les polygones traites sont convexes !
  915. DO 550 i = 1, 3
  916. vN1x = -vT1A(i,2)
  917. vN1y = vT1A(i,1)
  918. scal = vN1x * cPT1A(i,1) + vN1y * cPT1A(i,2)
  919. *
  920. ipos = 0
  921. ineg = 0
  922. jpos = 0
  923. jneg = 0
  924. do 551 j = 1, nPIn0
  925. test(j) = (vN1x * cPIn0(j,1) + vN1y * cPIn0(j,2)) - scal
  926. if (test(j).gt.1.d-10) then
  927. * if (test(j).gt.0.) then
  928. ipos = ipos + 1
  929. if (jneg.ne.0 .and. jpos.eq.0) jpos = j
  930. else if (test(j).lt.-1.d-10) then
  931. * else if (test(j).lt.0.) then
  932. ineg = ineg + 1
  933. if (jneg.eq.0) jneg = j
  934. else
  935. test(j) = 0.
  936. endif
  937. 551 continue
  938. *
  939. * 1) Cas ou ipos = 0 :
  940. * 1.1) il n'y pas d'intersection (ineg=nPIn0)
  941. * 1.2) l'intersection se limite a 1 point (ineg=nPIn0-1)
  942. * 1.3) l'intersection a 1 segment sur une frontiere (ineg=nPIn0-2)
  943. * La surface correspondante est nulle et il n'y a donc pas de contact !
  944. if (ipos.eq.0) then
  945. C*DBG write(ioimp,*) 'Intersection a ',nPIn0-ineg,' points'
  946. goto 520
  947. endif
  948. * 2) Cas ou ipos > 0 : Il y a une intersection a calculer
  949. * 2.1) ineg = 0 : Tous les points sont "du bon cote" et l'intersection
  950. * correspond au polygone de depart !
  951. if (ineg.eq.0) then
  952. C*DBG write(ioimp,*) 'Intersection du bon cote - rien a faire'
  953. goto 550
  954. endif
  955. * 2.2) ineg > 0 : Il y a deux intersections a determiner car il y a
  956. * au moins un point du "mauvais cote"
  957. jpos = max(1,jpos)
  958. jpr = jpos - 1
  959. if (jpos.eq.1) jpr = nPIn0
  960. nPIn = 1
  961. if (test(jpr).lt.0.) then
  962. r_z = test(jpos) / (test(jpos) - test(jpr))
  963. cPIn(nPIn,1) = cPIn0(jpos,1)
  964. & + r_z*(cPIn0(jpr,1) - cPIn0(jpos,1))
  965. cPIn(nPIn,2) = cPIn0(jpos,2)
  966. & + r_z*(cPIn0(jpr,2) - cPIn0(jpos,2))
  967. else
  968. cPIn(nPIn,1) = cPIn0(jpr,1)
  969. cPIn(nPIn,2) = cPIn0(jpr,2)
  970. endif
  971. do 552 j = 1, ipos
  972. nPIn = nPIn + 1
  973. cPIn(nPIn,1) = cPIn0(jpos,1)
  974. cPIn(nPIn,2) = cPIn0(jpos,2)
  975. jpr = jpos
  976. jpos = jpos + 1
  977. if (jpos.gt.nPIn0) jpos = 1
  978. 552 continue
  979. nPIn = nPIn + 1
  980. if (test(jpos).lt.0.) then
  981. r_z = test(jpr) / (test(jpr) - test(jpos))
  982. cPIn(nPIn,1) = cPIn0(jpr,1)
  983. & + r_z*(cPIn0(jpos,1) - cPIn0(jpr,1))
  984. cPIn(nPIn,2) = cPIn0(jpr,2)
  985. & + r_z*(cPIn0(jpos,2) - cPIn0(jpr,2))
  986. else
  987. cPIn(nPIn,1) = cPIn0(jpos,1)
  988. cPIn(nPIn,2) = cPIn0(jpos,2)
  989. endif
  990. * Mise a jour cPIn0 pour la suite du traitement des intersections...
  991. nPIn0 = nPIn
  992. do 554 j = 1, nPIn0
  993. cPIn0(j,1) = cPIn(j,1)
  994. cPIn0(j,2) = cPIn(j,2)
  995. 554 continue
  996. C*DBG-F write(ioimp,*) 'Intersection ',i,' a ',nPIn0,' points'
  997. C*DBG-F do j = 1, nPIn0
  998. C*DBG-F write(ioimp,*) ' ',j,nbpts0+6+j,cPIn0(j,1),cPIn0(j,2)
  999. C*DBG-F enddo
  1000. *
  1001. 550 CONTINUE
  1002. *
  1003. * Calcul de la surface de l'intersection et de son centre de gravite
  1004. r_z = cPIn0(nPIn0,1)*cPIn0(1,2) - cPIn0(1,1)*cPIn0(nPIn0,2)
  1005. SuPIn = r_z
  1006. xGIn = (cPIn0(nPIn0,1)+cPIn0(1,1)) * r_z
  1007. yGIn = (cPIn0(nPIn0,2)+cPIn0(1,2)) * r_z
  1008. do 560 i = 1, nPIn0-1
  1009. j = i + 1
  1010. r_z = cPIn0(i,1)*cPIn0(j,2) - cPIn0(j,1)*cPIn0(i,2)
  1011. SuPIn = SuPIn + r_z
  1012. xGIn = xGIn + (cPIn0(i,1)+cPIn0(j,1)) * r_z
  1013. yGIn = yGIn + (cPIn0(i,2)+cPIn0(j,2)) * r_z
  1014. 560 continue
  1015. if (SuPIn .lt. 1.E-7*max(SuT1A,SuT2A) ) then
  1016. C*DBG write(ioimp,*) 'Intersection a une surface negligeable !'
  1017. goto 520
  1018. endif
  1019. r_z = X1s3 / SupIn
  1020. SuPIn = 0.5 * SupIn
  1021. xGIn = xGIn * r_z
  1022. yGIn = yGIn * r_z
  1023. C*DBG write(ioimp,*) 'Surface Intersection :',SuPIn
  1024. *
  1025. * Calcul des coordonnees barycentriques du centre de gravite
  1026. * pour le triangle T1A
  1027. xpip1 = xGIn - cPT1A(1,1)
  1028. ypip1 = yGIn - cPT1A(1,2)
  1029. xpip2 = xGIn - cPT1A(2,1)
  1030. ypip2 = yGIn - cPT1A(2,2)
  1031. xpip3 = xGIn - cPT1A(3,1)
  1032. ypip3 = yGIn - cPT1A(3,2)
  1033. b1T1 = xpip2 * ypip3 - ypip2 * xpip3
  1034. b2T1 = xpip3 * ypip1 - ypip3 * xpip1
  1035. b3T1 = xpip1 * ypip2 - ypip1 * xpip2
  1036. bt = b1T1 + b2T1 + b3T1
  1037. if (abs(bt-SuT1A) .gt. 1.E-3)
  1038. & write(ioimp,*) 'Prob bt-SuT1A',bt,SuT1A
  1039. b1T1 = b1T1 / SuT1A
  1040. b2T1 = b2T1 / SuT1A
  1041. b3T1 = b3T1 / SuT1A
  1042. *
  1043. * pour le triangle T2A
  1044. xpip1 = xGIn - cPT2A(1,1)
  1045. ypip1 = yGIn - cPT2A(1,2)
  1046. xpip2 = xGIn - cPT2A(2,1)
  1047. ypip2 = yGIn - cPT2A(2,2)
  1048. xpip3 = xGIn - cPT2A(3,1)
  1049. ypip3 = yGIn - cPT2A(3,2)
  1050. b1T2 = xpip2 * ypip3 - ypip2 * xpip3
  1051. b2T2 = xpip3 * ypip1 - ypip3 * xpip1
  1052. b3T2 = xpip1 * ypip2 - ypip1 * xpip2
  1053. bt = b1T2 + b2T2 + b3T2
  1054. if (abs(bt-SuT2A) .gt. 1.E-3)
  1055. & write(ioimp,*) 'Prob bt-SuT2A',bt,SuT2A
  1056. b1T2 = b1T2 / SuT2A
  1057. b2T2 = b2T2 / SuT2A
  1058. b3T2 = b3T2 / SuT2A
  1059. *
  1060. * Calcul du jeu
  1061. xjeu = (b1T2*dq1C + b2T2*dq2C + b3T2*dq3C)
  1062. & - (b1T1*dp1C + b2T1*dp2C + b3T1*dp3C)
  1063. C*DBG write(ioimp,*) 'Jeu =',xjeu
  1064. *
  1065. * Ajout d'une relation triangle-triangle
  1066. *
  1067. nelri0 = nelri0 + 1
  1068. nelch = nelch + 1
  1069. *
  1070. * on ajuste les differents segments si necesaire
  1071. *
  1072. if (nelri0.gt.nelrig) then
  1073. nelrig = nelrig + incnel
  1074. segadj,xmatri
  1075. nbelem = nbelem + incnel
  1076. nbnn = 7
  1077. segadj,meleme
  1078. nbnn = 1
  1079. segadj,ipt7
  1080. n = n + incnel
  1081. segadj,mpoval
  1082. endif
  1083. *
  1084. * Mise a jour du meleme
  1085. *
  1086. if (ipt1.num(1,iel).eq.ILGNI) then
  1087. segact ipt1*mod
  1088. nbpts=xcoor(/1)/(idim+1)+1
  1089. segadj mcoord
  1090. ipt1.num(1,iel)=nbpts
  1091. if (ipt6.ne.0) ipt6.num(1,iel)=nbpts
  1092. endif
  1093. num(1,nelri0) = ipt1.num(1,iel)
  1094. num(2,nelri0) = ip1
  1095. num(3,nelri0) = ip2
  1096. num(4,nelri0) = ip3
  1097. num(5,nelri0) = jp1
  1098. num(6,nelri0) = jp2
  1099. num(7,nelri0) = jp3
  1100. icolor(nelri0) = 1
  1101. *$
  1102. * Mise a jour de xmatri
  1103. *
  1104. * lambda
  1105. re( 1,1,nelri0) = 0.d0
  1106. * ip1
  1107. re( 2,1,nelri0) = xnC * b1T1 * SuPIn
  1108. re( 3,1,nelri0) = ynC * b1T1 * SuPIn
  1109. re( 4,1,nelri0) = znC * b1T1 * SuPIn
  1110. * ip2
  1111. re( 5,1,nelri0) = xnC * b2T1 * SuPIn
  1112. re( 6,1,nelri0) = ynC * b2T1 * SuPin
  1113. re( 7,1,nelri0) = znC * b2T1 * SuPIn
  1114. * ip3
  1115. re( 8,1,nelri0) = xnC * b3T1 * SuPIn
  1116. re( 9,1,nelri0) = ynC * b3T1 * SuPIn
  1117. re(10,1,nelri0) = znC * b3T1 * SuPIn
  1118. * jp1
  1119. re(11,1,nelri0) = -xnC * b1T2 * SuPIn
  1120. re(12,1,nelri0) = -ynC * b1T2 * SuPIn
  1121. re(13,1,nelri0) = -znC * b1T2 * SuPIn
  1122. * jp2
  1123. re(14,1,nelri0) = -xnC * b2T2 * SuPIn
  1124. re(15,1,nelri0) = -ynC * b2T2 * SuPIn
  1125. re(16,1,nelri0) = -znC * b2T2 * SuPIn
  1126. * jp3
  1127. re(17,1,nelri0) = -xnC * b3T2 * SuPIn
  1128. re(18,1,nelri0) = -ynC * b3T2 * SuPIn
  1129. re(19,1,nelri0) = -znC * b3T2 * SuPIn
  1130. * on transpose
  1131. do 580 ic = 2, nligrp
  1132. re(1,ic,nelri0) = re(ic,1,nelri0)
  1133. 580 continue
  1134. * le reste est nul
  1135. *
  1136. * remplissage du champoint de depi (jeu)
  1137. *
  1138. ipt7.num(1,nelch) = ipt1.num(1,iel)
  1139. vpocha(nelch,1) = (xjeu - xjr) * SuPIn
  1140. C*DBG-F call prfaible(ipt1,iel,nPIn0,mfaible)
  1141. *
  1142. 520 CONTINUE
  1143. *
  1144. * Ajustement au plus juste puis desactivation des segments lies
  1145. * a la rigidite du type 0
  1146. if (nelri0.ne.nelrig) then
  1147. nelrig = nelri0
  1148. segadj,xmatri
  1149. nbelem = nelri0
  1150. nbnn = 7
  1151. segadj,meleme
  1152. endif
  1153. segdes,meleme,xmatri
  1154. *
  1155. c write(ioimp,*) ' nb relation type 0 ',nelri0
  1156. * if (nelri0.eq.0) irigel(6,nrigel)=0
  1157. *
  1158. * Destruction des segments locaux devenus inutiles
  1159. segsup,mfaible
  1160. *
  1161. C* GOTO 1000
  1162.  
  1163. *=======================================================================
  1164. *= Fin du sous-programme :
  1165. *=======================================================================
  1166. 1000 CONTINUE
  1167. *
  1168. * Ajustement au plus juste du chpoint de depi (jeu) : mpoval et ipt7
  1169. * puis desactivation du chpoint
  1170. if (vpocha(/1).ne.nelch) then
  1171. n = nelch
  1172. nc=1
  1173. segadj,mpoval
  1174. nbnn = 1
  1175. nbelem = nelch
  1176. nbsous = 0
  1177. nbref = 0
  1178. segadj,ipt7
  1179. endif
  1180. segdes,mpoval,ipt7
  1181. segdes,msoupo,mchpoi
  1182. * Desctivation de la matrice de raideur de contact
  1183. segdes,mrigid
  1184. *
  1185. * Reunion des relations portant sur le meme multiplicateur de lagrange
  1186. *
  1187. call impofu(MRIGID,MCHPOI)
  1188. *
  1189. * Ecriture des resultats de l'operateur :
  1190. call ecrobj('RIGIDITE',mrigid)
  1191. call ecrobj('CHPOINT ',mchpoi)
  1192. *
  1193. 900 CONTINUE
  1194. segdes,ipt1
  1195.  
  1196. return
  1197. end
  1198.  
  1199.  
  1200.  
  1201.  
  1202.  
  1203.  
  1204.  
  1205.  
  1206.  
  1207.  
  1208.  
  1209.  
  1210.  
  1211.  

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