Télécharger impo32.eso

Retour à la liste

Numérotation des lignes :

impo32
  1. C IMPO32 SOURCE PV090527 24/02/01 21:15:03 11829
  2.  
  3. * impo bloc en 3D
  4.  
  5. SUBROUTINE IMPO32(ipt6,ipt8,itcont)
  6.  
  7. IMPLICIT INTEGER(I-N)
  8. IMPLICIT REAL*8 (A-H,O-Z)
  9.  
  10.  
  11. -INC PPARAM
  12. -INC CCOPTIO
  13. -INC CCREEL
  14. -INC CCGEOME
  15.  
  16. -INC SMCOORD
  17. -INC SMELEME
  18. -INC SMRIGID
  19. -INC SMCHPOI
  20. -INC SMCHAML
  21.  
  22. * Contact en formulation forte :
  23. * Segments utiles dans les cas dits pathologiques
  24. segment mfopa1
  25. integer lpoin(nfopa1,3)
  26. real*8 xpoin(nfopa1,3)
  27. real*8 zje1(nfopa1)
  28. endsegment
  29. segment mfopa2
  30. integer lsegm(nfopa2,4)
  31. real*8 xsegm(nfopa2,3)
  32. real*8 zje2(nfopa2)
  33. endsegment
  34. * directions moyennes aux sommets
  35. segment mfopa3
  36. real*8 xms(nbpts),yms(nbpts),zms(nbpts)
  37. endsegment
  38. * positions compressees
  39. segment mfopa4
  40. real*8 tco(nbel1),xw(nbel1)
  41. integer ipnum(nbel1),ico(nbel1),iw(nbel1)
  42. endsegment
  43. segment mfopa5
  44. integer ind(indt)
  45. endsegment
  46.  
  47. * Contact en formulation faible :
  48. * Segment utile pour le calcul de l'intersection des 2 triangles
  49. * projetes dans le plan de contact intermediaire
  50. segment mfaible
  51. real*8 cPT1C(3,3),cPT2C(3,3), cPT1A(4,2),cPT2A(4,2)
  52. real*8 cPIn0(6,2), cPIn(6,2), test(6)
  53. real*8 vT1A(3,2), vT2A(3,2)
  54. real*8 SuT1A, SuT2A, SuPIn
  55. *DBG-F real*8 xGIn,yGIn, b1T1,b2T1,b3T1, b1T2,b2T2,b3T2
  56. endsegment
  57.  
  58. PARAMETER ( X1s3=0.333333333333333333333333333333333333333333D0 )
  59. PARAMETER ( sqr3=sqrt(3.d0))
  60.  
  61. character*4 modepl(3),moforc(3)
  62.  
  63. data modepl /'UX ','UY ','UZ '/
  64. data moforc /'FX ','FY ','FZ '/
  65. *
  66.  
  67. idimp1 = IDIM + 1
  68. *
  69. * Lecture du maillage support des conditions de contact
  70. * il s(agit la du premier maillage de contact, de type tri3
  71. *
  72. *
  73. call lirobj('MAILLAGE',ipt1,1,iretou)
  74. call actobj('MAILLAGE',ipt1,1)
  75. ** write(6,*) ' ipt1 dans impo32 ',ipt1
  76. ** write(6,*) ' ipt6 dans impo32 ',ipt6
  77. ** write(6,*) ' ipt8 dans impo32 ',ipt8
  78. isous=0
  79. call lirent(isous,0,iretou)
  80. mchel1=0
  81. call lirent(mchel1,0,iretou)
  82. C
  83. if (ierr.ne.0) return
  84. *
  85. * Activation du maillage et petites verifications
  86. *
  87. segact ipt8
  88. segact ipt6
  89. segact ipt1
  90. ** write(6,*) ' ipt1 2 dans impo32 ',ipt1
  91. nbno6 = ipt6.num(/1)
  92. nbel6 = ipt6.num(/2)
  93. nbno8 = ipt8.num(/1)
  94. nbel8 = ipt8.num(/2)
  95. nbel1 = ipt1.num(/2)
  96. ** write(6,*) 'nbel6 ',nbel6
  97. ** write(6,*) 'nbel1 ',nbel1
  98. if (ipt1.lisous(/1).ne.0) call erreur(25)
  99. if (ipt6.itypel.ne.4) call erreur(16)
  100. if (nbno6.ne.3 ) call erreur(16)
  101. if (ierr.ne.0) goto 900
  102. *
  103. * Increment sur le nombre d'elements de contact retenus et utilises
  104. * dans le chpoint (mpoval et igeoc) et rigidite pour agrandir les
  105. * longueurs des segments adequats
  106. incnel = 2000
  107. *
  108. * Creation de la raideur des conditions de contact
  109. * Remplissage de l'entete commun
  110. *
  111. nrigel = 1
  112. segini,mrigid
  113. ichole = 0
  114. imgeo1 = 0
  115. imgeo2 = 0
  116. isupeq = 0
  117. iforig = ifour
  118. coerig(1) = 1.
  119. mtymat='RIGIDITE'
  120. *
  121. * MCHAML materiau associe au MMODEL
  122. MELVA1 = 0
  123. MELVA2 = 0
  124. IF (MCHEL1.NE.0) THEN
  125. SEGACT, MCHEL1
  126. * recherche rapide d'un maillage correspondant dans le mchaml
  127. DO 210 n2 = 1,mchel1.imache(/1)
  128. * write(ioimp,*) ' n2 imache(n2) ipt1 ',n2,mchel1.imache(n2),ipt1
  129. if (mchel1.imache(n2).eq.ipt1) goto 220
  130. 210 continue
  131.  
  132. goto 230
  133. 220 continue
  134. MCHAM1 = MCHEL1.ICHAML(n2)
  135. SEGACT, MCHAM1
  136. NBCO = MCHAM1.NOMCHE(/2)
  137. DO JCMP = 1,NBCO
  138. IF (MCHAM1.NOMCHE(JCMP).EQ.'JEU') THEN
  139. MELVA1 = MCHAM1.IELVAL(JCMP)
  140. SEGACT, MELVA1
  141. NELJ = MELVA1.VELCHE(/2)
  142. NPTELJ = min(MELVA1.VELCHE(/1),4)
  143. C
  144. C Utilise pour le zonage
  145. xjmax = -xpetit
  146. xjmin = xpetit
  147. DO IJI2=1,NELJ
  148. DO IJI1=1,MELVA1.VELCHE(/1)
  149. xjmax=max(xjmax,MELVA1.VELCHE(IJI1,IJI2))
  150. xjmin=min(xjmin,MELVA1.VELCHE(IJI1,IJI2))
  151. ENDDO
  152. ENDDO
  153. C
  154. ENDIF
  155. IF (MCHAM1.NOMCHE(JCMP).EQ.'ADHE') THEN
  156. MELVA2 = MCHAM1.IELVAL(JCMP)
  157. SEGACT, MELVA2
  158. NELA = MELVA2.VELCHE(/2)
  159. NPTELA = min(MELVA2.VELCHE(/1),4)
  160. ENDIF
  161. ENDDO
  162. ENDIF
  163. 230 continue
  164. *
  165. * Creation du chpoint de depi
  166. *
  167. nat=1
  168. nsoupo=1
  169. segini mchpoi
  170. mtypoi='DEPIMP'
  171. mochde='engendré par impose'
  172. ifopoi=ifour
  173. jattri(1)=2
  174. *
  175. nc=2
  176. IF (MELVA2.NE.0) nc=3
  177. segini msoupo
  178. ipchp(1)=msoupo
  179. nocomp(1)='FLX '
  180. nocomp(2)='SCAL'
  181. IF (MELVA2.NE.0) THEN
  182. nocomp(3)='FADH'
  183. ENDIF
  184. *
  185. nbnn =1
  186. nbelem=incnel
  187. nbref =0
  188. nbsous=0
  189. segini ipt7
  190. ipt7.itypel=1
  191. igeoc=ipt7
  192. *
  193. n=incnel
  194. segini mpoval
  195. ipoval=mpoval
  196. *
  197. * Nombre de noeuds dans le chpoint (en totalite) : ipt7 et mpoval
  198. nelch = 0
  199.  
  200. *=======================================================================
  201. * Formulation "forte" (standard) du contact :
  202. *=======================================================================
  203. * Element de contact 3D a 3 noeuds
  204. ** write(6,*) ' itcont dans impo32 ',itcont
  205. * itcont 2: formulation faible
  206. if (itcont.eq.2) goto 500
  207.  
  208. * Calcul des directions moyennes de la boite d'encadrement, et de la taille max
  209. segact mcoord
  210. segini mfopa3
  211.  
  212. xmin=xgrand
  213. xmax=-xgrand
  214. ymin=xgrand
  215. ymax=-xgrand
  216. zmin=xgrand
  217. zmax=-xgrand
  218. tamax=0.d0
  219.  
  220. do 820 iel6=1,nbel6
  221. *
  222. ip1=ipt6.num(1,iel6)
  223. ip2=ipt6.num(2,iel6)
  224. ip3=ipt6.num(3,iel6)
  225. ipv1 = (ip1-1)*idimp1
  226. xp1 = xcoor(ipv1+1)
  227. yp1 = xcoor(ipv1+2)
  228. zp1 = xcoor(ipv1+3)
  229. ipv2 = (ip2-1)*idimp1
  230. xp2 = xcoor(ipv2+1)
  231. yp2 = xcoor(ipv2+2)
  232. zp2 = xcoor(ipv2+3)
  233. ipv3 = (ip3-1)*idimp1
  234. xp3 = xcoor(ipv3+1)
  235. yp3 = xcoor(ipv3+2)
  236. zp3 = xcoor(ipv3+3)
  237.  
  238. xmin=min(xmin,xp1,xp2,xp3)
  239. xmax=max(xmax,xp1,xp2,xp3)
  240. ymin=min(ymin,yp1,yp2,yp3)
  241. ymax=max(ymax,yp1,yp2,yp3)
  242. zmin=min(zmin,zp1,zp2,zp3)
  243. zmax=max(zmax,zp1,zp2,zp3)
  244.  
  245. *
  246. * normale au plan (123)
  247. *
  248. x12 = xp2 - xp1
  249. y12 = yp2 - yp1
  250. z12 = zp2 - zp1
  251. x23 = xp3 - xp2
  252. y23 = yp3 - yp2
  253. z23 = zp3 - zp2
  254. xn = y12*z23 - z12*y23
  255. yn = z12*x23 - x12*z23
  256. zn = x12*y23 - y12*x23
  257. sn = xn*xn + yn*yn + zn*zn
  258. sn = max(sqrt(abs(sn)),xpetit*1d10)
  259. tamax=max(sn,tamax)
  260.  
  261. xn = xn/sn
  262. yn = yn/sn
  263. zn = zn/sn
  264. xms(ip1)=xms(ip1)+xn
  265. yms(ip1)=yms(ip1)+yn
  266. zms(ip1)=zms(ip1)+zn
  267. xms(ip2)=xms(ip2)+xn
  268. yms(ip2)=yms(ip2)+yn
  269. zms(ip2)=zms(ip2)+zn
  270. xms(ip3)=xms(ip3)+xn
  271. yms(ip3)=yms(ip3)+yn
  272. zms(ip3)=zms(ip3)+zn
  273. 820 continue
  274. tamax = tamax * 5.
  275.  
  276. do 822 iel8=1,nbel8
  277. *
  278. ip1=ipt8.num(1,iel8)
  279. ip2=ipt8.num(2,iel8)
  280. ip3=ipt8.num(3,iel8)
  281. ipv1 = (ip1-1)*idimp1
  282. xp1 = xcoor(ipv1+1)
  283. yp1 = xcoor(ipv1+2)
  284. zp1 = xcoor(ipv1+3)
  285. ipv2 = (ip2-1)*idimp1
  286. xp2 = xcoor(ipv2+1)
  287. yp2 = xcoor(ipv2+2)
  288. zp2 = xcoor(ipv2+3)
  289. ipv3 = (ip3-1)*idimp1
  290. xp3 = xcoor(ipv3+1)
  291. yp3 = xcoor(ipv3+2)
  292. zp3 = xcoor(ipv3+3)
  293. *
  294. * normale au plan (123)
  295. *
  296. x12 = xp2 - xp1
  297. y12 = yp2 - yp1
  298. z12 = zp2 - zp1
  299. x23 = xp3 - xp2
  300. y23 = yp3 - yp2
  301. z23 = zp3 - zp2
  302. xn = y12*z23 - z12*y23
  303. yn = z12*x23 - x12*z23
  304. zn = x12*y23 - y12*x23
  305. sn = xn*xn + yn*yn + zn*zn
  306. sn = max(sqrt(abs(sn)),xpetit*1d10)
  307. xn = xn/sn
  308. yn = yn/sn
  309. zn = zn/sn
  310. xms(ip1)=xms(ip1)+xn
  311. yms(ip1)=yms(ip1)+yn
  312. zms(ip1)=zms(ip1)+zn
  313. xms(ip2)=xms(ip2)+xn
  314. yms(ip2)=yms(ip2)+yn
  315. zms(ip2)=zms(ip2)+zn
  316. xms(ip3)=xms(ip3)+xn
  317. yms(ip3)=yms(ip3)+yn
  318. zms(ip3)=zms(ip3)+zn
  319. 822 continue
  320. do 821 ip = 1,nbpts
  321. sn = xms(ip)*xms(ip)+yms(ip)*yms(ip)+zms(ip)*zms(ip)
  322. sn = max(sqrt(abs(sn)),xpetit*1d10)
  323. xms(ip)=xms(ip)/sn
  324. yms(ip)=yms(ip)/sn
  325. zms(ip)=zms(ip)/sn
  326. ** write(6,*) 'ip ',ip,xms(ip),yms(ip),zms(ip)
  327. 821 continue
  328. if(ipt1.itypel.ne.22) call erreur(16)
  329. ** write(6,*) 'ipt8.itypel ',ipt8.itypel, ipt7
  330. if(ierr.ne.0) return
  331. *
  332. * Nombre d'iterations pour la detection de la direction du contact
  333. nbiter=16
  334. *
  335. * Relation du type 3 : noeud-triangle
  336. *-------------------------------------
  337. *
  338. * creation du meleme associe a la relation
  339. *
  340. nbsous=0
  341. nbref =0
  342. nbnn =5
  343. nbelem=incnel
  344. segini meleme
  345. itypel=22
  346. irigel(1,nrigel)=meleme
  347. *
  348. * creation du descriptif commun a toutes les raideurs
  349. *
  350. nligrp = 13
  351. nligrd = nligrp
  352. segini,descr
  353. lisinc(1)='LX '
  354. lisdua(1)='FLX '
  355. noelep(1)=1
  356. noeled(1)=1
  357. do 10 i=2, nligrp, 3
  358. lisinc(i) =modepl(1)
  359. lisinc(i+1)=modepl(2)
  360. lisinc(i+2)=modepl(3)
  361. lisdua(i) =moforc(1)
  362. lisdua(i+1)=moforc(2)
  363. lisdua(i+2)=moforc(3)
  364. noelep(i) =(i+4)/3
  365. noelep(i+1)=noelep(i)
  366. noelep(i+2)=noelep(i)
  367. noeled(i) =noelep(i)
  368. noeled(i+1)=noelep(i)
  369. noeled(i+2)=noelep(i)
  370. 10 continue
  371. segdes,descr
  372. irigel(3,nrigel) = descr
  373. *
  374. * creation du segment xmatri
  375. *
  376. nelrig = incnel
  377. segini,xmatri
  378. irigel(4,nrigel) = xmatri
  379. *
  380. * ce qu'on cree est unilateral
  381. *
  382. irigel(6,nrigel)=1
  383. *
  384. * ce qu'on cree est symetrique
  385. *
  386. irigel(7,nrigel)=0
  387. *
  388. * Nombre d'elements dans la rigidite du type 3
  389. nelri3 = 0
  390. *
  391. * Preparation du zonage de ipt1
  392. *
  393. nbzx=(xmax-xmin)/tamax
  394. nbzy=(ymax-ymin)/tamax
  395. nbzz=(zmax-zmin)/tamax
  396. *
  397. * preconditionnement de ipt1: mfopa4
  398. *
  399. segini mfopa4
  400. xmult=3.1415926*(xmax-xmin)
  401. ymult=2.7182818*(ymax-ymin)
  402. zmult=1.*(zmax-zmin)
  403. tmult=sqrt(xmult**2+ymult**2+zmult**2)
  404. xmult=xmult/tmult
  405. ymult=ymult/tmult
  406. zmult=zmult/tmult
  407. tmin=xgrand
  408. tmax=-xgrand
  409. do j=1,nbel1
  410. ip=ipt1.num(2,j)
  411. ipv = (ip-1)*idimp1
  412. xp = xcoor(ipv+1)
  413. yp = xcoor(ipv+2)
  414. zp = xcoor(ipv+3)
  415. tco(j)=xp*xmult+yp*ymult+zp*zmult
  416. tmin=min(tco(j),tmin)
  417. tmax=max(tco(j),tmax)
  418. ipnum(j)=ip
  419. ico(j)=j
  420. enddo
  421. nbzt=(tmax-tmin)/tamax
  422. nbzt=max(nbzt,1)
  423. nbzt=min(nbel1,nbzt)
  424. ** write(6,*) ' nbzt ',nbzt
  425.  
  426. *
  427. * trier selon x
  428. *
  429. call triflo(tco,xw,ico,iw,nbel1)
  430. *
  431. * indexer
  432. *
  433. indt=nbzt+1
  434. segini mfopa5
  435. do i=nbel1,1,-1
  436. id=nbzt*(tco(i)-tmin)/(tmax-tmin)+1
  437. ind(id)=i
  438. enddo
  439. do i=1,nbzt
  440. if (ind(i+1).eq.0) ind(i+1)=ind(i)
  441. if(ind(i+1).lt.ind(i)) call erreur(5)
  442. enddo
  443.  
  444.  
  445. * Boucle sur le maillage support du contact(frottement)
  446. *
  447. ** write(6,*) ' boucle 19 nbel6 ',nbel6
  448. DO 19 iel6=1,nbel6
  449. *
  450. xjr = 0d0
  451. ip1=ipt6.num(1,iel6)
  452. ip2=ipt6.num(2,iel6)
  453. ip3=ipt6.num(3,iel6)
  454. * Recuperation des coordonees des noeuds du triangle
  455. ipv = (ip1-1)*idimp1
  456. xp1 = xcoor(ipv+1)
  457. yp1 = xcoor(ipv+2)
  458. zp1 = xcoor(ipv+3)
  459. ipv = (ip2-1)*idimp1
  460. xp2 = xcoor(ipv+1)
  461. yp2 = xcoor(ipv+2)
  462. zp2 = xcoor(ipv+3)
  463. ipv = (ip3-1)*idimp1
  464. xp3 = xcoor(ipv+1)
  465. yp3 = xcoor(ipv+2)
  466. zp3 = xcoor(ipv+3)
  467. * Centre de gravite du triangle
  468. xg = (xp1+xp2+xp3) /3.d0
  469. yg = (yp1+yp2+yp3) /3.d0
  470. zg = (zp1+zp2+zp3) /3.d0
  471. * critere de distance
  472. d1 = (xg-xp1)**2 + (yg-yp1)**2 + (zg-zp1)**2
  473. d2 = (xg-xp2)**2 + (yg-yp2)**2 + (zg-zp2)**2
  474. d3 = (xg-xp3)**2 + (yg-yp3)**2 + (zg-zp3)**2
  475. dist2 = max(d1,d2,d3)
  476. sqdist = sqrt(dist2)
  477. sqdist4 = sqdist * 4
  478. * Triangle un peu plus grand pour les tests
  479. scale=1.00 + xszpre
  480. xp1e=xg+(xp1-xg)*scale
  481. yp1e=yg+(yp1-yg)*scale
  482. zp1e=zg+(zp1-zg)*scale
  483. xp2e=xg+(xp2-xg)*scale
  484. yp2e=yg+(yp2-yg)*scale
  485. zp2e=zg+(zp2-zg)*scale
  486. xp3e=xg+(xp3-xg)*scale
  487. yp3e=yg+(yp3-yg)*scale
  488. zp3e=zg+(zp3-zg)*scale
  489. 25 continue
  490. * rechercher le noeud a tester dans le deuxieme maillage qui est sous forme mult avec en
  491. * deuxieme position le point physique
  492. ** write(6,*) ' boucle 20 ipt1 ',ipt1.num(/2)
  493. xgm = xg -sqdist4
  494. xgp = xg +sqdist4
  495. ygm = yg -sqdist4
  496. ygp = yg +sqdist4
  497. zgm = zg -sqdist4
  498. zgp = zg +sqdist4
  499. *
  500. * calcul zone du centre de gravite
  501. *
  502. tc=xg*xmult+yg*ymult+zg*zmult
  503. *
  504. xmn = (xms(ip1)+xms(ip2)+xms(ip3))*X1s3
  505. ymn = (yms(ip1)+yms(ip2)+yms(ip3))*X1s3
  506. zmn = (zms(ip1)+zms(ip2)+zms(ip3))*X1s3
  507. xzonag = sqdist4*sqr3
  508. tjeup=0.d0
  509. tjeum=0.d0
  510.  
  511. if (MELVA1.ne.0) then
  512. tjeup=xmult*(xjmax*xmn)+ymult*(xjmax*ymn)+zmult*(xjmax*zmn)
  513. tjeum=xmult*(xjmin*xmn)+ymult*(xjmin*ymn)+zmult*(xjmin*zmn)
  514. endif
  515. *
  516. izg=nbzt*((tc-xzonag-tjeup)-tmin)/(tmax-tmin)+1
  517. izg=max(izg,1)
  518. izg=min(izg,indt)
  519. * debut de zone
  520. indb=ind(izg)
  521. do 20 iz=indb,nbel1
  522. if(tco(iz).lt.(tc-xzonag-tjeup)) goto 20
  523. if(tco(iz).gt.(tc+xzonag-tjeum)) then
  524. * write(6,*) ' sortie pour ',iz, ' en ',iz-indb+1
  525. goto 18
  526. endif
  527. iel1 = ico(iz)
  528. jp = ipnum(iel1)
  529. * Verification que pas relation du point sur lui meme
  530. if (jp.eq.ip1) goto 20
  531. if (jp.eq.ip2) goto 20
  532. if (jp.eq.ip3) goto 20
  533. * verification rapide en norme L1 d'abord
  534. ipv = (jp-1)*idimp1
  535. xp = xcoor(ipv+1)
  536. yp = xcoor(ipv+2)
  537. zp = xcoor(ipv+3)
  538. *
  539. xpp=xp
  540. ypp=yp
  541. zpp=zp
  542. if (MELVA1.ne.0) then
  543. nel1 = min (iel1,nelj)
  544. xjr = melva1.velche(nptelj,nel1)
  545. xpp=xpp-xjr*xmn
  546. ypp=ypp-xjr*ymn
  547. zpp=zpp-xjr*zmn
  548. endif
  549. if(xpp.lt.xgm.or.xpp.gt.xgp) goto 20
  550. if(ypp.lt.ygm.or.ypp.gt.ygp) goto 20
  551. if(zpp.lt.zgm.or.zpp.gt.zgp) goto 20
  552. *
  553. * Verification si autre point dans la zone de selection
  554. * verif distance au centre de gravite
  555. dp = ((xg-xpp)**2 + (yg-ypp)**2 + (zg-zpp)**2)
  556. ** write(6,*) 'dp dist2',dp,dist2,xjr
  557. *** if (dp.gt.xzonag**2) then
  558. *** goto 20
  559. *** endif
  560. C*DBG write(ioimp,*) 'contact test distance ok',dp,d1,d2,d3
  561. dist = SQRT(dist2)
  562.  
  563. * verif position par rapport aux cotes
  564.  
  565. * cote 1 2
  566. x12 = xp2 - xp1
  567. y12 = yp2 - yp1
  568. z12 = zp2 - zp1
  569. sv = sqrt(x12**2+y12**2+z12**2)
  570. xv=x12/sv
  571. yv=y12/sv
  572. zv=z12/sv
  573. * normale locale (1 et 2)
  574. xnl=xms(ip1)+xms(ip2)
  575. ynl=yms(ip1)+yms(ip2)
  576. znl=zms(ip1)+zms(ip2)
  577.  
  578. * vecteur reference
  579. xn=y12*znl-z12*ynl
  580. yn=z12*xnl-x12*znl
  581. zn=x12*ynl-y12*xnl
  582. dn = sqrt(xn*xn+yn*yn+zn*zn)
  583. scal = (xpp-xp1e)*xv + (ypp-yp1e)*yv + (zpp-zp1e)*zv
  584. xm = xp1e + scal*xv
  585. ym = yp1e + scal*yv
  586. zm = zp1e + scal*zv
  587. dpm = sqrt(abs((xpp-xm)**2+(ypp-ym)**2+(zpp-zm)**2))
  588. scal = (xpp-xm)*xn + (ypp-ym)*yn + (zpp-zm)*zn
  589. ** if (dpm.lt.dist*xszpre) write(6,*) ' pt sur 1 2',dpm
  590. if (dpm.gt.dist*xszpre.and.scal.gt.0.707d0*dn*dpm) then
  591. * write(6,*) ' 1 dpm scal ',dpm,scal
  592. goto 20
  593. endif
  594. dpm=max(xpetit,dpm)
  595. * write(ioimp,*) 'contact test position 1 ok',
  596. * & scal/(dn*dpm),scal,dn,dpm
  597. *
  598. * cote 2 3
  599. x23 = xp3 - xp2
  600. y23 = yp3 - yp2
  601. z23 = zp3 - zp2
  602. sv = sqrt(x23**2+y23**2+z23**2)
  603. xv=x23/sv
  604. yv=y23/sv
  605. zv=z23/sv
  606. * normale locale (2 et 3)
  607. xnl=xms(ip2)+xms(ip3)
  608. ynl=yms(ip2)+yms(ip3)
  609. znl=zms(ip2)+zms(ip3)
  610.  
  611. * vecteur reference
  612. xn=y23*znl-z23*ynl
  613. yn=z23*xnl-x23*znl
  614. zn=x23*ynl-y23*xnl
  615. dn = sqrt(xn*xn+yn*yn+zn*zn)
  616. scal = (xpp-xp2e)*xv + (ypp-yp2e)*yv + (zpp-zp2e)*zv
  617. xm = xp2e + scal*xv
  618. ym = yp2e + scal*yv
  619. zm = zp2e + scal*zv
  620. dpm = sqrt(abs((xpp-xm)**2+(ypp-ym)**2+(zpp-zm)**2))
  621. scal = (xpp-xm)*xn + (ypp-ym)*yn + (zpp-zm)*zn
  622. ** if (dpm.lt.dist*xszpre) write(6,*) ' pt sur 2 3',dpm
  623. if (dpm.gt.dist*xszpre.and.scal.gt.0.707d0*dn*dpm) then
  624. * write(6,*) ' 2 dpm scal ',dpm,scal
  625. goto 20
  626. endif
  627. dpm=max(xpetit,dpm)
  628. * write(ioimp,*) 'contact test position 2 ok',
  629. * & scal/(dn*dpm),scal,dn,dpm
  630. *
  631. * cote 3 1
  632. x31 = xp1 - xp3
  633. y31 = yp1 - yp3
  634. z31 = zp1 - zp3
  635. sv = sqrt(x31**2+y31**2+z31**2)
  636. xv=x31/sv
  637. yv=y31/sv
  638. zv=z31/sv
  639. * normale locale (3 et 1)
  640. xnl=xms(ip3)+xms(ip1)
  641. ynl=yms(ip3)+yms(ip1)
  642. znl=zms(ip3)+zms(ip1)
  643. * vecteur reference
  644. xn=y31*znl-z31*ynl
  645. yn=z31*xnl-x31*znl
  646. zn=x31*ynl-y31*xnl
  647. dn = sqrt(xn*xn+yn*yn+zn*zn)
  648. scal = (xpp-xp3e)*xv + (ypp-yp3e)*yv + (zpp-zp3e)*zv
  649. xm = xp3e + scal*xv
  650. ym = yp3e + scal*yv
  651. zm = zp3e + scal*zv
  652. dpm = sqrt(abs((xpp-xm)**2+(ypp-ym)**2+(zpp-zm)**2))
  653. scal = (xpp-xm)*xn + (ypp-ym)*yn + (zpp-zm)*zn
  654. ** if (dpm.lt.dist*xszpre) write(6,*) ' pt sur 3 1',dpm
  655. if (dpm.gt.dist*xszpre.and.scal.gt.0.707d0*dn*dpm) then
  656. * write(6,*) ' 3 dpm scal ',dpm,scal
  657. goto 20
  658. endif
  659. dpm=max(xpetit,dpm)
  660. * write(ioimp,*) 'contact test position 2 ok',
  661. * & scal/(dn*dpm),scal,dn,dpm
  662. *
  663. * on a un point ou imposer la relation
  664. C*DBG write(ioimp,*) 'contact potentiel '
  665. * Quelle est la relation ???
  666. *
  667. * direction de la relation = normale au plan (123)
  668. *
  669. * normale reelle
  670. xnr = y12*z23 - z12*y23
  671. ynr = z12*x23 - x12*z23
  672. znr = x12*y23 - y12*x23
  673. dnr = sqrt(xnr*xnr+ynr*ynr+znr*znr)
  674. if (dnr.lt.dist*xszpre.AND.iimpi.ne.0) write(ioimp,*) ' pb dnr'
  675. xnr = xnr / dnr
  676. ynr = ynr / dnr
  677. znr = znr / dnr
  678. * normale ponderee
  679. xn=xms(ip1)+xms(ip2)+xms(ip3)
  680. yn=yms(ip1)+yms(ip2)+yms(ip3)
  681. zn=zms(ip1)+zms(ip2)+zms(ip3)
  682. dn = sqrt(xn*xn+yn*yn+zn*zn)
  683. if (.not.(dn.lt.1d-10).and. .not.(dn.gt.1d-10).AND.iimpi.ne.0)
  684. & write(ioimp,*) 'Prob 4.1 - impo32'
  685. if (abs(dn).le.xpetit.AND.iimpi.ne.0)
  686. & write(ioimp,*) 'Prob 4.2 - impo32'
  687. xn = xn / dn
  688. yn = yn / dn
  689. zn = zn / dn
  690. *
  691. * Distance (jeu) du point jp au plan de contact
  692. * Puis calcul de la projection sur le plan de contact
  693. *
  694. ** write(6,*) 'ip1 xms ',ip1,xms(ip1),yms(ip1),zms(ip1)
  695. ** write(6,*) 'ip2 xms ',ip2,xms(ip2),yms(ip2),zms(ip2)
  696. ** write(6,*) 'ip3 xms ',ip3,xms(ip3),yms(ip3),zms(ip3)
  697. ** write(6,*) 'xn ',xn,yn,zn
  698. iter = 0
  699. xjeu = (xp-xg)*xnr + (yp-yg)*ynr + (zp-zg)*znr
  700. 21 continue
  701. angn = xn * xnr + yn*ynr + zn*znr
  702. if (angn.le.xpetit.AND.iimpi.ne.-1)
  703. & write(ioimp,*) 'angn negatif ',angn
  704. if(angn.le.xpetit) goto 20
  705. xjeuv = xjeu / angn
  706. *
  707. * Recherche de ses coordonnées barycentriques
  708. * qui sont les surfaces des sous triangles
  709. *
  710. xq = xp - xjeuv*xn
  711. yq = yp - xjeuv*yn
  712. zq = zp - xjeuv*zn
  713. xb1 = (yq-yp2)*(zq-zp3)-(zq-zp2)*(yq-yp3)
  714. yb1 = (zq-zp2)*(xq-xp3)-(xq-xp2)*(zq-zp3)
  715. zb1 = (xq-xp2)*(yq-yp3)-(yq-yp2)*(xq-xp3)
  716. xb2 = (yq-yp3)*(zq-zp1)-(zq-zp3)*(yq-yp1)
  717. yb2 = (zq-zp3)*(xq-xp1)-(xq-xp3)*(zq-zp1)
  718. zb2 = (xq-xp3)*(yq-yp1)-(yq-yp3)*(xq-xp1)
  719. xb3 = (yq-yp1)*(zq-zp2)-(zq-zp1)*(yq-yp2)
  720. yb3 = (zq-zp1)*(xq-xp2)-(xq-xp1)*(zq-zp2)
  721. zb3 = (xq-xp1)*(yq-yp2)-(yq-yp1)*(xq-xp2)
  722. b1 = xb1*xnr + yb1*ynr + zb1*znr
  723. b2 = xb2*xnr + yb2*ynr + zb2*znr
  724. b3 = xb3*xnr + yb3*ynr + zb3*znr
  725. bt = b1 + b2 + b3
  726. * normalement pas utile
  727. * element retourne a cause des grands deplacements
  728. if (bt.lt.xpetit) then
  729. write (ioimp,*) ' bt negatif dans impo32 '
  730. call soucis(719)
  731. ** bt=xpetit*2.d0
  732. goto 20
  733. endif
  734. if (bt.lt.0.d0.AND.iimpi.ne.0) then
  735. write (ioimp,*) ' bt negatif dans impo32 '
  736. print *,'xp1 yp1 zp1',xp1,yp1,zp1
  737. print *,'xp2 yp2 zp2',xp2,yp2,zp2
  738. print *,'xp3 yp3 zp3',xp3,yp3,zp3
  739. print *,'xp yp zp',xp,yp,zp
  740. print *,'xn yn zn',xn,yn,zn
  741. print *,'b1 b2 b3 bt',b1,b2,b3,bt
  742. * goto 20
  743. endif
  744. bsurf = bt / 2
  745. if (.not.(bt.lt.1d-10) .and. .not.(bt.gt.1d-10).AND.iimpi.ne.0)
  746. & write(ioimp,*) 'Prob 5.1 - impo32'
  747. if (abs(bt).le.xpetit) write(ioimp,*) 'Prob 5.2 - impo32'
  748. b1 = b1 / bt
  749. b2 = b2 / bt
  750. b3 = b3 / bt
  751. * recalcul de la direction et du jeu en fonction des normales aux sommets
  752. * on arete l'interpolation de la normale au bord de l'element
  753. bb1=max(b1,0.d0)
  754. bb2=max(b2,0.d0)
  755. bb3=max(b3,0.d0)
  756. bm = bb1 + bb2 + bb3
  757. bb1 = bb1 / bm
  758. bb2 = bb2 / bm
  759. bb3 = bb3 / bm
  760. xnp = xn
  761. ynp = yn
  762. znp = zn
  763. xn = bb1*xms(ip1)+ bb2*xms(ip2)+ bb3*xms(ip3)
  764. yn = bb1*yms(ip1)+ bb2*yms(ip2)+ bb3*yms(ip3)
  765. zn = bb1*zms(ip1)+ bb2*zms(ip2)+ bb3*zms(ip3)
  766. sn = sqrt (xn*xn + yn*yn + zn*zn)
  767. xn=xn/sn
  768. yn=yn/sn
  769. zn=zn/sn
  770. diff=abs(xn-xnp)+abs(yn-ynp)+abs(zn-znp)
  771. * write (6,*) ' iter b* *n ',iter,b1,b2,b3,xn,yn,zn
  772. * recalcul en fonction de la nouvelle normale
  773. iter=iter+1
  774. if (iter.gt.64) then
  775. ** write(6,*) ' impo32 diff ',diff
  776. goto 20
  777. endif
  778. if (diff.gt.1d-10) goto 21
  779. ** write(6,*) ' impo32 iter ',iter
  780. ** write (6,*) ' iter b* *n ',iter,b1,b2,b3,xn,yn,zn
  781.  
  782.  
  783. * si on a deja traverse, les trois coordonnees barycentriques doivent etre positives
  784. ** if (xjeuv-xjr.lt.-xszpre*dist) then
  785. ** if (b1.lt.-xszpre) goto 20
  786. ** if (b2.lt.-xszpre) goto 20
  787. ** if (b3.lt.-xszpre) goto 20
  788. ** endif
  789. * Si on est en dehors, on projette sur l'arete (ou pas)
  790. ** if (b1.lt.0.d0.or.b2.lt.0.d0.or.b3.lt.0.d0) then
  791. ** xq=xp1*bb1+xp2*bb2+xp3*bb3
  792. ** yq=yp1*bb1+yp2*bb2+yp3*bb3
  793. ** zq=zp1*bb1+zp2*bb2+zp3*bb3
  794. ** xnn=xp-xq
  795. ** ynn=yp-yq
  796. ** znn=zp-zq
  797. ** snn=sqrt(abs(xnn*xnn+ynn*ynn+znn*znn))
  798. ** if (snn.lt.xszpre*dist) write (6,*) ' snn petit ',snn
  799. * sinon on prend la direction reelle
  800. ** if (snn.gt.xszpre*dist) then
  801. ** xn=xnn/sn
  802. ** yn=ynn/sn
  803. ** zn=znn/sn
  804. ** endif
  805. ** endif
  806. xjeu1 = (xp-xp1)*xn + (yp-yp1)*yn + (zp-zp1)*zn
  807. xjeu2 = (xp-xp2)*xn + (yp-yp2)*yn + (zp-zp2)*zn
  808. xjeu3 = (xp-xp3)*xn + (yp-yp3)*yn + (zp-zp3)*zn
  809. *****pv xjeuv = xjeu1 * bb1 + xjeu2 * bb2 + xjeu3 * bb3
  810. ** write (6,*) 'b1 b2 b3',b1,b2,b3
  811. ** write (6,*) 'xjeu xjeuv ',xjeu,xjeuv,xjeu1,xjeu2,xjeu3
  812. xjeu = xjeuv - xjr
  813. ** write (6,*) ' normale ',xn,yn,zn,' jeu ',xjeu1,xjeu2,xjeu3
  814. * verif bon cote (a la tolerance pres)
  815. ** write (6,*) ' xjeu ',xjeu,dist,iel6
  816. ** if (xjeu.lt.-1*dist) goto 20
  817. ** if (xjeu.gt.3*dist) goto 20
  818. * verif compatibilite avec la normale au poin impactant
  819. if (xms(jp)*xn+yms(jp)*yn+zms(jp)*zn.gt. -0.0d0) then
  820. * write (6,*) ' impo32 normales incompatibes 1 ',
  821. * > jp,xjeu,dpm,dist
  822. goto 20
  823. endif
  824.  
  825.  
  826.  
  827. C*DBG write(ioimp,*) ' b1 b2 b3 ',b1,b2,b3
  828. 1954 continue
  829. * points a l'exterieur?
  830. * on met une ponderation et un rayon d'acceptation
  831. pond = (1D4+1) - 1D4 * bm
  832. if (pond.le.0.d0) goto 20
  833. pond=max(pond,0.d0)
  834. pond=min(pond,1.d0)
  835.  
  836. ** if (xjeu.lt.0.d0)
  837. ** > write (6,*) 'pt traverse',jp,b1,b2,b3,xjeuv-xjr,dist
  838.  
  839. *
  840. * Ajout d'une relation noeud-triangle
  841. *
  842. nelri3 = nelri3 + 1
  843. nelch = nelch + 1
  844. *
  845. * on ajuste les differents segments si necesaire
  846. *
  847. if (nelri3.gt.nelrig) then
  848. nelrig = nelrig + incnel
  849. segadj,xmatri
  850. nbelem = nbelem + incnel
  851. nbnn = 5
  852. segadj,meleme
  853. nbnn = 1
  854. segadj,ipt7
  855. n = n + incnel
  856. segadj,mpoval
  857. endif
  858. *
  859. * Mise a jour du meleme
  860. ** write(6,*) 'ipt8.num(/2)',ipt8.num(/2)
  861. num(1,nelri3) = ipt1.num(1,iel1)
  862. num(2,nelri3) = ip1
  863. num(3,nelri3) = ip2
  864. num(4,nelri3) = ip3
  865. num(5,nelri3) = jp
  866. ** write(6,*) ' nouveau meleme',num(1,nelri3),num(2,nelri3),
  867. ** > num(3,nelri3),num(4,nelri3),num(5,nelri3) *
  868. * Mise a jour de xmatri
  869. * lambda
  870. re( 1,1,nelri3) = 0d0
  871. * ip1
  872. re( 2,1,nelri3) = xn * bb1 * bsurf * pond
  873. re( 3,1,nelri3) = yn * bb1 * bsurf * pond
  874. re( 4,1,nelri3) = zn * bb1 * bsurf * pond
  875. * ip2
  876. re( 5,1,nelri3) = xn * bb2 * bsurf * pond
  877. re( 6,1,nelri3) = yn * bb2 * bsurf * pond
  878. re( 7,1,nelri3) = zn * bb2 * bsurf * pond
  879. * ip3
  880. re( 8,1,nelri3) = xn * bb3 * bsurf * pond
  881. re( 9,1,nelri3) = yn * bb3 * bsurf * pond
  882. re(10,1,nelri3) = zn * bb3 * bsurf * pond
  883. * jp
  884. re(11,1,nelri3) = -xn * bsurf * pond
  885. re(12,1,nelri3) = -yn * bsurf * pond
  886. re(13,1,nelri3) = -zn * bsurf * pond
  887. * on transpose
  888. do 30 ic = 2, nligrp
  889. re(1,ic,nelri3) = re(ic,1,nelri3)
  890. 30 continue
  891. * le reste est nul
  892. *
  893. * remplissage du champoint de depi (jeu)
  894. *
  895. ipt7.num(1,nelch) = ipt1.num(1,iel1)
  896. ** write(6,*) ' ipt1.num(1,iel1) ',ipt1.num(1,iel1)
  897. vpocha(nelch,1) = xjeu * bsurf * pond
  898. vpocha(nelch,2) = bsurf * pond
  899. IF (MELVA2.NE.0) THEN
  900. NEL2 = min (iel1,NELA)
  901. VPOCHA(nelch,3) = max(MELVA2.VELCHE(NPTELA,NEL2),0.D0)*
  902. > bsurf * pond
  903. ENDIF
  904. ** write(ioimp,*) ' jeu ',xjeu,ip1,ip2,ip3,jp,b1,b2,b3,xn,yn,zn
  905. *
  906. 20 CONTINUE
  907. 18 CONTINUE
  908. 19 CONTINUE
  909. *
  910. * Ajustement au plus juste puis desactivation des segments lies
  911. * a la rigidite du type 3
  912. if (nelri3.ne.nelrig) then
  913. nelrig = nelri3
  914. segadj,xmatri
  915. nbelem = nelri3
  916. nbnn = 5
  917. segadj,meleme
  918. endif
  919. segdes,xmatri
  920. *
  921. ** write(ioimp,*) ' nb relation type 3 ',nelri3
  922. ** write(6,fmt='(10i8)') (num(5,nelr),nelr=1,nelri3)
  923.  
  924.  
  925.  
  926.  
  927.  
  928. C*DBG if (nelri3.eq.0) irigel(6,nrigel)=0
  929. *
  930. * Fin du traitement de la formulation standard/forte
  931. * ----------------------------------------------------
  932. 300 CONTINUE
  933. * Destruction des segments locaux devenus inutiles
  934. segsup mfopa3,mfopa4,mfopa5
  935.  
  936. GOTO 1000
  937.  
  938. *=======================================================================
  939. *= Formulation "faible" 3D : relation triangle-triangle
  940. *=======================================================================
  941. * Element de contact 3D a 7 noeuds (1+6)
  942. 500 CONTINUE
  943. *
  944. * Petit segment de travail
  945. segini,mfaible
  946. *
  947. * Relation du type 0 : triangle-triangle (faible)
  948. *---------------------------------------
  949. *
  950. * Creation du meleme associe a la relation
  951. *
  952. nbnn = 7
  953. nbelem = incnel
  954. nbsous = 0
  955. nbref = 0
  956. segini,meleme
  957. itypel = 22
  958. irigel(1,nrigel) = meleme
  959. *
  960. * Creation du descriptif commun a toutes les raideurs
  961. *
  962. nligrp = 19
  963. nligrd = nligrp
  964. segini,descr
  965. lisinc(1) = 'LX '
  966. lisdua(1) = 'FLX '
  967. noelep(1) = 1
  968. noeled(1) = 1
  969. do 510 i = 2, nligrp, 3
  970. lisinc(i ) = modepl(1)
  971. lisinc(i+1) = modepl(2)
  972. lisinc(i+2) = modepl(3)
  973. lisdua(i ) = moforc(1)
  974. lisdua(i+1) = moforc(2)
  975. lisdua(i+2) = moforc(3)
  976. noelep(i ) = (i+5)/3
  977. noelep(i+1) = noelep(i)
  978. noelep(i+2) = noelep(i)
  979. noeled(i ) = noelep(i)
  980. noeled(i+1) = noelep(i)
  981. noeled(i+2) = noelep(i)
  982. 510 continue
  983. segdes,descr
  984. irigel(3,nrigel) = descr
  985. *
  986. * creation du segment xmatri
  987. *
  988. nelrig = incnel
  989. segini,xmatri
  990. irigel(4,nrigel) = xmatri
  991. *
  992. * ce qu'on cree est unilateral
  993. *
  994. irigel(6,nrigel) = 1
  995. *
  996. * ce qu'on cree est symetrique
  997. *
  998. irigel(7,nrigel) = 0
  999. *
  1000. * Nombre d'elements crees dans meleme=irigel(1,nrigel), ipt7 et mpoval
  1001. nelri0 = 0
  1002. *
  1003. * Boucle sur les elements du maillage de contact/frottement "faible"
  1004. *
  1005. DO 519 iel1 = 1, nbel1
  1006. *
  1007. xjr = 0d0
  1008. ip1 = ipt1.num(2,iel1)
  1009. ip2 = ipt1.num(3,iel1)
  1010. ip3 = ipt1.num(4,iel1)
  1011. * Definition du triangle T1(ip1,ip2,ip3)
  1012. * - Coordonnees des noeuds de T1
  1013. ipv = (ip1-1)*idimp1
  1014. xp1 = xcoor(ipv+1)
  1015. yp1 = xcoor(ipv+2)
  1016. zp1 = xcoor(ipv+3)
  1017. ipv = (ip2-1)*idimp1
  1018. xp2 = xcoor(ipv+1)
  1019. yp2 = xcoor(ipv+2)
  1020. zp2 = xcoor(ipv+3)
  1021. ipv = (ip3-1)*idimp1
  1022. xp3 = xcoor(ipv+1)
  1023. yp3 = xcoor(ipv+2)
  1024. zp3 = xcoor(ipv+3)
  1025. * - Barycentre de T1
  1026. xgT1 = (xp1 + xp2 + xp3) * X1s3
  1027. ygT1 = (yp1 + yp2 + yp3) * X1s3
  1028. zgT1 = (zp1 + zp2 + zp3) * X1s3
  1029. * - Normale au triangle T1
  1030. xnT1 = (yp1-yp2)*(zp2-zp3) - (zp1-zp2)*(yp2-yp3)
  1031. ynT1 = (zp1-zp2)*(xp2-xp3) - (xp1-xp2)*(zp2-zp3)
  1032. znT1 = (xp1-xp2)*(yp2-yp3) - (yp1-yp2)*(xp2-xp3)
  1033. dnT1 = SQRT(xnT1*xnT1+ynT1*ynT1+znT1*znT1)
  1034. IF (.not.(dnT1.lt.1d-10) .and. .not.(dnT1.gt.1d-10).AND.
  1035. & iimpi.ne.0) write(ioimp,*) 'FAIBle - 1.1 - impo32'
  1036. IF (abs(dnT1).le.xpetit) write(ioimp,*) 'FAIBle - 1.2 - impo32'
  1037. xnT1 = xnT1 / dnT1
  1038. ynT1 = ynT1 / dnT1
  1039. znT1 = znT1 / dnT1
  1040. C*DBG write(ioimp,*) 'Normale au plan T1',xnT1,ynT1,znT1,dnT1
  1041. * - Rayon de la sphere "englobante" de l'element (centre=barycentre)
  1042. d1 = (xgT1-xp1)**2 + (ygT1-yp1)**2 + (zgT1-zp1)**2
  1043. d2 = (xgT1-xp2)**2 + (ygT1-yp2)**2 + (zgT1-zp2)**2
  1044. d3 = (xgT1-xp3)**2 + (ygT1-yp3)**2 + (zgT1-zp3)**2
  1045. RayT1 = SQRT(max(d1,d2,d3))
  1046. DO 520 iel6 = 1, nbel6
  1047. if (MELVA1.ne.0) then
  1048. nel1 = min (iel6,nelj)
  1049. xjr = melva1.velche(nptelj,nel1)
  1050. endif
  1051. jp1 = ipt6.num(1,iel6)
  1052. * - Les noeuds 2 et 3 sont intervertis (fait aupravant dans impo31)
  1053. * Cela permet aux triangles T1 et T2 d etre orientes dans le meme sens.
  1054. jp3 = ipt6.num(2,iel6)
  1055. jp2 = ipt6.num(3,iel6)
  1056. C*DBG write(ioimp,*) iel,ip1,ip2,ip3,jp1,jp2,jp3,ipt1.num(1,iel)
  1057. *
  1058. * Verification (provisoire) qu'il n'y a pas de noeud double
  1059. if (ip1.eq.jp1) goto 520
  1060. if (ip1.eq.jp2) goto 520
  1061. if (ip1.eq.jp3) goto 520
  1062. if (ip2.eq.jp1) goto 520
  1063. if (ip2.eq.jp2) goto 520
  1064. if (ip2.eq.jp3) goto 520
  1065. if (ip3.eq.jp1) goto 520
  1066. if (ip3.eq.jp2) goto 520
  1067. if (ip3.eq.jp3) goto 520
  1068. *
  1069. *
  1070. * Definition du triangle T2(jp1,jp2,jp3)
  1071. * - Coordonnees des noeuds de T2
  1072. ipv = (jp1-1)*idimp1
  1073. xq1 = xcoor(ipv+1)
  1074. yq1 = xcoor(ipv+2)
  1075. zq1 = xcoor(ipv+3)
  1076. ipv = (jp2-1)*idimp1
  1077. xq2 = xcoor(ipv+1)
  1078. yq2 = xcoor(ipv+2)
  1079. zq2 = xcoor(ipv+3)
  1080. ipv = (jp3-1)*idimp1
  1081. xq3 = xcoor(ipv+1)
  1082. yq3 = xcoor(ipv+2)
  1083. zq3 = xcoor(ipv+3)
  1084. * - Barycentre de T2
  1085. xgT2 = (xq1 + xq2 + xq3) * X1s3
  1086. ygT2 = (yq1 + yq2 + yq3) * X1s3
  1087. zgT2 = (zq1 + zq2 + zq3) * X1s3
  1088. * - Normale au triangle T2
  1089. xnT2 = (yq1-yq2)*(zq2-zq3) - (zq1-zq2)*(yq2-yq3)
  1090. ynT2 = (zq1-zq2)*(xq2-xq3) - (xq1-xq2)*(zq2-zq3)
  1091. znT2 = (xq1-xq2)*(yq2-yq3) - (yq1-yq2)*(xq2-xq3)
  1092. dnT2 = SQRT(xnT2*xnT2+ynT2*ynT2+znT2*znT2)
  1093. IF (.not.(dnT2.lt.1d-10) .and. .not.(dnT2.gt.1d-10).AND.
  1094. & iimpi.ne.0) write(ioimp,*) 'FAIBle - 2.1 - impo32'
  1095. IF (abs(dnT2).le.xpetit.AND.iimpi.ne.0)
  1096. & write(ioimp,*) 'FAIBle - 2.2 - impo32'
  1097. xnT2 = xnT2 / dnT2
  1098. ynT2 = ynT2 / dnT2
  1099. znT2 = znT2 / dnT2
  1100. C*DBG write(ioimp,*) 'Normale au plan T2',xnT2,ynT2,znT2,dnT2
  1101. * - Rayon de la sphere "englobante" de l'element (centre=barycentre)
  1102. d1 = (xgT2-xq1)**2 + (ygT2-yq1)**2 + (zgT2-zq1)**2
  1103. d2 = (xgT2-xq2)**2 + (ygT2-yq2)**2 + (zgT2-zq2)**2
  1104. d3 = (xgT2-xq3)**2 + (ygT2-yq3)**2 + (zgT2-zq3)**2
  1105. RayT2 = SQRT(max(d1,d2,d3))
  1106.  
  1107. * Orientation respective correcte des 2 normales
  1108. * Les triangles T1 et T2 etant decrits dans le meme sens ("prisme"),
  1109. * le produit scalaire de leur normale doit etre positif ou nul !
  1110. scal = xnT1 * xnT2 + ynT1 * ynT2 + znT1 * znT2
  1111. C*DBG write(ioimp,*) 'Prod.scal des normales',scal
  1112. if (scal.lt.0.) goto 520
  1113. *
  1114. * Distance entre les centres de gravite
  1115. * Les spheres "englobantes" doivent s'intersecter pour contact potentiel
  1116. dist = SQRT((xgT2-xgT1)**2 + (ygT2-ygT1)**2 + (zgT2-zgT1)**2)
  1117. * IF (dist*4.GT.RayT1+RayT2) goto 520
  1118. *
  1119. * Definition du plan de contact (surface intermediaire)
  1120. * - Normale (unitaire) du plan de contact
  1121. xnC = xnT1 + xnT2
  1122. ynC = ynT1 + ynT2
  1123. znC = znT1 + znT2
  1124. dnC = SQRT(xnC*xnC + ynC*ynC + znC*znC)
  1125. IF (.not.(dnC.lt.1d-10).and. .not.(dnC.GT.1D-10).AND.iimpi.ne.0)
  1126. & write(ioimp,*) 'FAIBle - 3.1 - impo32'
  1127. IF (abs(dnC).le.xpetit.AND.iimpi.ne.0)
  1128. & write(ioimp,*) 'FAIBle - 3.2 - impo32'
  1129. xnC = xnC / dnC
  1130. ynC = ynC / dnC
  1131. znC = znC / dnC
  1132. C*DBG write(ioimp,*) 'Normale au plan Contact',xnC,ynC,znC,dnC
  1133. *
  1134. * Criteres d'acceptation de l'element de contact
  1135. RayTT = RayT1 + RayT2
  1136. * - Les barycentres des triangles T1 et T2 sont-ils suffisamment proches :
  1137. * distance suivant la direction de contact (prise en compte du jeu)
  1138. VXB12 = xgT2 - xgT1
  1139. VYB12 = ygT2 - ygT1
  1140. VZB12 = zgT2 - zgT1
  1141. distn = (VXB12*xnC)+(VYB12*ynC)+(VZB12*znC)
  1142. test1 = distn
  1143. if (MELVA1.NE.0) then
  1144. test1 = test1 - xjr
  1145. endif
  1146. if (test1.GT.RayTT) GOTO 520
  1147. *
  1148. * distance dans un plan dont la normale est la direction de contact
  1149. distt = SQRT(dist**2 - distn**2)
  1150. if (distt.GT.RayTT) GOTO 520
  1151.  
  1152. * - Point definissant le plan de contact
  1153. * = milieu des barycentres des triangles T1 et T2
  1154. xgC = (xgT1 + xgT2) * 0.5
  1155. ygC = (ygT1 + ygT2) * 0.5
  1156. zgC = (zgT1 + zgT2) * 0.5
  1157. * - "Distance" de reference au plan de contact
  1158. dgC = xgC*xnC + ygC*ynC + zgC*znC
  1159. *
  1160. * Distance "signee" des triangles T1 et T2 au plan de contact
  1161. dp1C = (xp1*xnC + yp1*ynC + zp1*znC) - dgC
  1162. dp2C = (xp2*xnC + yp2*ynC + zp2*znC) - dgC
  1163. dp3C = (xp3*xnC + yp3*ynC + zp3*znC) - dgC
  1164. dq1C = (xq1*xnC + yq1*ynC + zq1*znC) - dgC
  1165. dq2C = (xq2*xnC + yq2*ynC + zq2*znC) - dgC
  1166. dq3C = (xq3*xnC + yq3*ynC + zq3*znC) - dgC
  1167. * Projection des triangles T1 et T2 sur le plan de contact
  1168. * T1C = triangle T1 projete sur plan de Contact
  1169. cPT1C(1,1) = xp1 - dp1C * xnC
  1170. cPT1C(1,2) = yp1 - dp1C * ynC
  1171. cPT1C(1,3) = zp1 - dp1C * znC
  1172. cPT1C(2,1) = xp2 - dp2C * xnC
  1173. cPT1C(2,2) = yp2 - dp2C * ynC
  1174. cPT1C(2,3) = zp2 - dp2C * znC
  1175. cPT1C(3,1) = xp3 - dp3C * xnC
  1176. cPT1C(3,2) = yp3 - dp3C * ynC
  1177. cPT1C(3,3) = zp3 - dp3C * znC
  1178. * T2C = triangle T2 projete sur plan de Contact
  1179. cPT2C(1,1) = xq1 - dq1C * xnC
  1180. cPT2C(1,2) = yq1 - dq1C * ynC
  1181. cPT2C(1,3) = zq1 - dq1C * znC
  1182. cPT2C(2,1) = xq2 - dq2C * xnC
  1183. cPT2C(2,2) = yq2 - dq2C * ynC
  1184. cPT2C(2,3) = zq2 - dq2C * znC
  1185. cPT2C(3,1) = xq3 - dq3C * xnC
  1186. cPT2C(3,2) = yq3 - dq3C * ynC
  1187. cPT2C(3,3) = zq3 - dq3C * znC
  1188. *
  1189. * On determine quelle est la composante maximale de la normale pour
  1190. * projeter les triangles T1C et T2C sur un plan parallele aux axes
  1191. * ("le plus orthogonal" a cette normale), ce qui maximise leur aire.
  1192. r_x = abs(xnC)
  1193. r_y = abs(ynC)
  1194. r_z = abs(znC)
  1195. if (r_x.gt.r_y) then
  1196. inC = 1
  1197. if (r_x.lt.r_z) inC = 3
  1198. else
  1199. inC = 2
  1200. if (r_y.lt.r_z) inC = 3
  1201. endif
  1202. * Projection des triangles T1C et T2C sur ce plan // aux axes (A)
  1203. * T1A = triangle T1C projete sur ce plan (A)
  1204. * T2A = triangle T2C projete sur ce plan (A)
  1205. * On prend soin d'orienter les triangles T1A et T2A dans le sens
  1206. * direct (en tenant compte du signe de la composante normale !)
  1207. GOTO (531,532,533), inC
  1208. call erreur(5)
  1209. return
  1210. * Normale NC selon "x" : A =(y,z)
  1211. 531 if (xnC.lt.0.) inC = -inC
  1212. inX = 2
  1213. inY = 3
  1214. goto 534
  1215. * Normale NC selon "y" : A =(x,z)
  1216. 532 if (ynC.lt.0.) inC = -inC
  1217. inX = 1
  1218. inY = 3
  1219. goto 534
  1220. * Normale NC selon "z" : A =(x,y)
  1221. 533 if (znC.lt.0.) inC = -inC
  1222. inX = 1
  1223. inY = 2
  1224. goto 534
  1225. 534 continue
  1226. cPT1A(1,1) = cPT1C(1,inX)
  1227. cPT1A(1,2) = cPT1C(1,inY)
  1228. cPT2A(1,1) = cPT2C(1,inX)
  1229. cPT2A(1,2) = cPT2C(1,inY)
  1230. if (inC.gt.0) then
  1231. cPT1A(2,1) = cPT1C(2,inX)
  1232. cPT1A(2,2) = cPT1C(2,inY)
  1233. cPT2A(2,1) = cPT2C(2,inX)
  1234. cPT2A(2,2) = cPT2C(2,inY)
  1235. cPT1A(3,1) = cPT1C(3,inX)
  1236. cPT1A(3,2) = cPT1C(3,inY)
  1237. cPT2A(3,1) = cPT2C(3,inX)
  1238. cPT2A(3,2) = cPT2C(3,inY)
  1239. else
  1240. cPT1A(2,1) = cPT1C(3,inX)
  1241. cPT1A(2,2) = cPT1C(3,inY)
  1242. cPT2A(2,1) = cPT2C(3,inX)
  1243. cPT2A(2,2) = cPT2C(3,inY)
  1244. cPT1A(3,1) = cPT1C(2,inX)
  1245. cPT1A(3,2) = cPT1C(2,inY)
  1246. cPT2A(3,1) = cPT2C(2,inX)
  1247. cPT2A(3,2) = cPT2C(2,inY)
  1248. endif
  1249. cPT1A(4,1) = cPT1A(1,1)
  1250. cPT1A(4,2) = cPT1A(1,2)
  1251. * Calcul des cotes des triangles T1A et T2A (12,23,31)
  1252. * Utile pour connaitre ensuite la normale aux cotes des triangles
  1253. vT1A(1,1) = cPT1A(2,1) - cPT1A(1,1)
  1254. vT1A(1,2) = cPT1A(2,2) - cPT1A(1,2)
  1255. vT1A(2,1) = cPT1A(3,1) - cPT1A(2,1)
  1256. vT1A(2,2) = cPT1A(3,2) - cPT1A(2,2)
  1257. vT1A(3,1) = cPT1A(1,1) - cPT1A(3,1)
  1258. vT1A(3,2) = cPT1A(1,2) - cPT1A(3,2)
  1259. C*nu vT2A(1,1) = cPT2A(2,1) - cPT2A(1,1)
  1260. vT2A(1,2) = cPT2A(2,2) - cPT2A(1,2)
  1261. C*nu vT2A(2,1) = cPT2A(3,1) - cPT2A(2,1)
  1262. vT2A(2,2) = cPT2A(3,2) - cPT2A(2,2)
  1263. C*nu vT2A(3,1) = cPT2A(1,1) - cPT2A(3,1)
  1264. vT2A(3,2) = cPT2A(1,2) - cPT2A(3,2)
  1265. * Calcul de la surface de chacun des triangles T1A et T2A
  1266. * (en fait, on calcule le double, mais par la suite on ne considere que
  1267. * des rapports de surfaces de triangle)
  1268. SuT1A = (cPT1A(2,1)+cPT1A(1,1)) * vT1A(1,2)
  1269. & + (cPT1A(3,1)+cPT1A(2,1)) * vT1A(2,2)
  1270. & + (cPT1A(1,1)+cPT1A(3,1)) * vT1A(3,2)
  1271. SuT2A = (cPT2A(2,1)+cPT2A(1,1)) * vT2A(1,2)
  1272. & + (cPT2A(3,1)+cPT2A(2,1)) * vT2A(2,2)
  1273. & + (cPT2A(1,1)+cPT2A(3,1)) * vT2A(3,2)
  1274. C*DBG write(ioimp,*) 'Surfaces T1A - T2A :', 0.5*SuT1A,0.5*SuT2A
  1275. C*DBG if (SuT1A.gt.100.*SuT2A .or. SuT2A.gt.100.*SuT1A)
  1276. C*DBG & write(ioimp,*) 'Rapport des surfaces tres important !'
  1277.  
  1278. * On initialise l'intersection des triangles avec T2A
  1279. nPIn0 = 3
  1280. do 540 i = 1, nPIn0
  1281. cPIn0(i,1) = cPT2A(i,1)
  1282. cPIn0(i,2) = cPT2A(i,2)
  1283. 540 continue
  1284. * Determination de l'intersection des triangles T1A et T2A
  1285. * en regardant progressivement l'intersection avec les cotes de T1A
  1286. * Note : On utilise le fait que les polygones traites sont convexes !
  1287. DO 550 i = 1, 3
  1288. vN1x = -vT1A(i,2)
  1289. vN1y = vT1A(i,1)
  1290. scal = vN1x * cPT1A(i,1) + vN1y * cPT1A(i,2)
  1291. *
  1292. ipos = 0
  1293. ineg = 0
  1294. jpos = 0
  1295. jneg = 0
  1296. do 551 j = 1, nPIn0
  1297. test(j) = (vN1x * cPIn0(j,1) + vN1y * cPIn0(j,2)) - scal
  1298. if (test(j).gt.1.d-10) then
  1299. * if (test(j).gt.0.) then
  1300. ipos = ipos + 1
  1301. if (jneg.ne.0 .and. jpos.eq.0) jpos = j
  1302. else if (test(j).lt.-1.d-10) then
  1303. * else if (test(j).lt.0.) then
  1304. ineg = ineg + 1
  1305. if (jneg.eq.0) jneg = j
  1306. else
  1307. test(j) = 0.
  1308. endif
  1309. 551 continue
  1310. *
  1311. * 1) Cas ou ipos = 0 :
  1312. * 1.1) il n'y pas d'intersection (ineg=nPIn0)
  1313. * 1.2) l'intersection se limite a 1 point (ineg=nPIn0-1)
  1314. * 1.3) l'intersection a 1 segment sur une frontiere (ineg=nPIn0-2)
  1315. * La surface correspondante est nulle et il n'y a donc pas de contact !
  1316. if (ipos.eq.0) then
  1317. C*DBG write(ioimp,*) 'Intersection a ',nPIn0-ineg,' points'
  1318. goto 520
  1319. endif
  1320. * 2) Cas ou ipos > 0 : Il y a une intersection a calculer
  1321. * 2.1) ineg = 0 : Tous les points sont "du bon cote" et l'intersection
  1322. * correspond au polygone de depart !
  1323. if (ineg.eq.0) then
  1324. C*DBG write(ioimp,*) 'Intersection du bon cote - rien a faire'
  1325. goto 550
  1326. endif
  1327. * 2.2) ineg > 0 : Il y a deux intersections a determiner car il y a
  1328. * au moins un point du "mauvais cote"
  1329. jpos = max(1,jpos)
  1330. jpr = jpos - 1
  1331. if (jpos.eq.1) jpr = nPIn0
  1332. nPIn = 1
  1333. if (test(jpr).lt.0.) then
  1334. r_z = test(jpos) / (test(jpos) - test(jpr))
  1335. cPIn(nPIn,1) = cPIn0(jpos,1)
  1336. & + r_z*(cPIn0(jpr,1) - cPIn0(jpos,1))
  1337. cPIn(nPIn,2) = cPIn0(jpos,2)
  1338. & + r_z*(cPIn0(jpr,2) - cPIn0(jpos,2))
  1339. else
  1340. cPIn(nPIn,1) = cPIn0(jpr,1)
  1341. cPIn(nPIn,2) = cPIn0(jpr,2)
  1342. endif
  1343. do 552 j = 1, ipos
  1344. nPIn = nPIn + 1
  1345. cPIn(nPIn,1) = cPIn0(jpos,1)
  1346. cPIn(nPIn,2) = cPIn0(jpos,2)
  1347. jpr = jpos
  1348. jpos = jpos + 1
  1349. if (jpos.gt.nPIn0) jpos = 1
  1350. 552 continue
  1351. nPIn = nPIn + 1
  1352. if (test(jpos).lt.0.) then
  1353. r_z = test(jpr) / (test(jpr) - test(jpos))
  1354. cPIn(nPIn,1) = cPIn0(jpr,1)
  1355. & + r_z*(cPIn0(jpos,1) - cPIn0(jpr,1))
  1356. cPIn(nPIn,2) = cPIn0(jpr,2)
  1357. & + r_z*(cPIn0(jpos,2) - cPIn0(jpr,2))
  1358. else
  1359. cPIn(nPIn,1) = cPIn0(jpos,1)
  1360. cPIn(nPIn,2) = cPIn0(jpos,2)
  1361. endif
  1362. * Mise a jour cPIn0 pour la suite du traitement des intersections...
  1363. nPIn0 = nPIn
  1364. do 554 j = 1, nPIn0
  1365. cPIn0(j,1) = cPIn(j,1)
  1366. cPIn0(j,2) = cPIn(j,2)
  1367. 554 continue
  1368. C*DBG-F write(ioimp,*) 'Intersection ',i,' a ',nPIn0,' points'
  1369. C*DBG-F do j = 1, nPIn0
  1370. C*DBG-F write(ioimp,*) ' ',j,nbpts0+6+j,cPIn0(j,1),cPIn0(j,2)
  1371. C*DBG-F enddo
  1372. *
  1373. 550 CONTINUE
  1374. *
  1375. * Calcul de la surface de l'intersection et de son centre de gravite
  1376. r_z = cPIn0(nPIn0,1)*cPIn0(1,2) - cPIn0(1,1)*cPIn0(nPIn0,2)
  1377. SuPIn = r_z
  1378. xGIn = (cPIn0(nPIn0,1)+cPIn0(1,1)) * r_z
  1379. yGIn = (cPIn0(nPIn0,2)+cPIn0(1,2)) * r_z
  1380. do 560 i = 1, nPIn0-1
  1381. j = i + 1
  1382. r_z = cPIn0(i,1)*cPIn0(j,2) - cPIn0(j,1)*cPIn0(i,2)
  1383. SuPIn = SuPIn + r_z
  1384. xGIn = xGIn + (cPIn0(i,1)+cPIn0(j,1)) * r_z
  1385. yGIn = yGIn + (cPIn0(i,2)+cPIn0(j,2)) * r_z
  1386. 560 continue
  1387. if (SuPIn .lt. 1.E-7*max(SuT1A,SuT2A) ) then
  1388. C*DBG write(ioimp,*) 'Intersection a une surface negligeable !'
  1389. goto 520
  1390. endif
  1391. r_z = X1s3 / SupIn
  1392. SuPIn = 0.5 * SupIn
  1393. xGIn = xGIn * r_z
  1394. yGIn = yGIn * r_z
  1395. C*DBG write(ioimp,*) 'Surface Intersection :',SuPIn
  1396. *
  1397. * Calcul des coordonnees barycentriques du centre de gravite
  1398. * pour le triangle T1A
  1399. xpip1 = xGIn - cPT1A(1,1)
  1400. ypip1 = yGIn - cPT1A(1,2)
  1401. xpip2 = xGIn - cPT1A(2,1)
  1402. ypip2 = yGIn - cPT1A(2,2)
  1403. xpip3 = xGIn - cPT1A(3,1)
  1404. ypip3 = yGIn - cPT1A(3,2)
  1405. b1T1 = xpip2 * ypip3 - ypip2 * xpip3
  1406. b2T1 = xpip3 * ypip1 - ypip3 * xpip1
  1407. b3T1 = xpip1 * ypip2 - ypip1 * xpip2
  1408. bt = b1T1 + b2T1 + b3T1
  1409. if (abs(bt-SuT1A) .gt. 1.E-3)
  1410. & write(ioimp,*) 'Prob bt-SuT1A',bt,SuT1A
  1411. b1T1 = b1T1 / SuT1A
  1412. b2T1 = b2T1 / SuT1A
  1413. b3T1 = b3T1 / SuT1A
  1414. *
  1415. * pour le triangle T2A
  1416. xpip1 = xGIn - cPT2A(1,1)
  1417. ypip1 = yGIn - cPT2A(1,2)
  1418. xpip2 = xGIn - cPT2A(2,1)
  1419. ypip2 = yGIn - cPT2A(2,2)
  1420. xpip3 = xGIn - cPT2A(3,1)
  1421. ypip3 = yGIn - cPT2A(3,2)
  1422. b1T2 = xpip2 * ypip3 - ypip2 * xpip3
  1423. b2T2 = xpip3 * ypip1 - ypip3 * xpip1
  1424. b3T2 = xpip1 * ypip2 - ypip1 * xpip2
  1425. bt = b1T2 + b2T2 + b3T2
  1426. if (abs(bt-SuT2A) .gt. 1.E-3)
  1427. & write(ioimp,*) 'Prob bt-SuT2A',bt,SuT2A
  1428. b1T2 = b1T2 / SuT2A
  1429. b2T2 = b2T2 / SuT2A
  1430. b3T2 = b3T2 / SuT2A
  1431. *
  1432. * Calcul du jeu
  1433. xjeu = (b1T2*dq1C + b2T2*dq2C + b3T2*dq3C)
  1434. & - (b1T1*dp1C + b2T1*dp2C + b3T1*dp3C)
  1435. C*DBG write(ioimp,*) 'Jeu =',xjeu
  1436. *
  1437. * Ajout d'une relation triangle-triangle
  1438. *
  1439. nelri0 = nelri0 + 1
  1440. nelch = nelch + 1
  1441. *
  1442. * on ajuste les differents segments si necesaire
  1443. *
  1444. if (nelri0.gt.nelrig) then
  1445. nelrig = nelrig + incnel
  1446. segadj,xmatri
  1447. nbelem = nbelem + incnel
  1448. nbnn = 7
  1449. segadj,meleme
  1450. nbnn = 1
  1451. segadj,ipt7
  1452. n = n + incnel
  1453. segadj,mpoval
  1454. endif
  1455. *
  1456. * Mise a jour du meleme
  1457. *
  1458. num(1,nelri0) = ipt1.num(1,iel1)
  1459. num(2,nelri0) = ip1
  1460. num(3,nelri0) = ip2
  1461. num(4,nelri0) = ip3
  1462. num(5,nelri0) = jp1
  1463. num(6,nelri0) = jp2
  1464. num(7,nelri0) = jp3
  1465. icolor(nelri0) = 1
  1466. *$
  1467. * Mise a jour de xmatri
  1468. *
  1469. * lambda
  1470. re( 1,1,nelri0) = 0.d0
  1471. * ip1
  1472. re( 2,1,nelri0) = xnC * b1T1 * SuPIn
  1473. re( 3,1,nelri0) = ynC * b1T1 * SuPIn
  1474. re( 4,1,nelri0) = znC * b1T1 * SuPIn
  1475. * ip2
  1476. re( 5,1,nelri0) = xnC * b2T1 * SuPIn
  1477. re( 6,1,nelri0) = ynC * b2T1 * SuPin
  1478. re( 7,1,nelri0) = znC * b2T1 * SuPIn
  1479. * ip3
  1480. re( 8,1,nelri0) = xnC * b3T1 * SuPIn
  1481. re( 9,1,nelri0) = ynC * b3T1 * SuPIn
  1482. re(10,1,nelri0) = znC * b3T1 * SuPIn
  1483. * jp1
  1484. re(11,1,nelri0) = -xnC * b1T2 * SuPIn
  1485. re(12,1,nelri0) = -ynC * b1T2 * SuPIn
  1486. re(13,1,nelri0) = -znC * b1T2 * SuPIn
  1487. * jp2
  1488. re(14,1,nelri0) = -xnC * b2T2 * SuPIn
  1489. re(15,1,nelri0) = -ynC * b2T2 * SuPIn
  1490. re(16,1,nelri0) = -znC * b2T2 * SuPIn
  1491. * jp3
  1492. re(17,1,nelri0) = -xnC * b3T2 * SuPIn
  1493. re(18,1,nelri0) = -ynC * b3T2 * SuPIn
  1494. re(19,1,nelri0) = -znC * b3T2 * SuPIn
  1495. * on transpose
  1496. do 580 ic = 2, nligrp
  1497. re(1,ic,nelri0) = re(ic,1,nelri0)
  1498. 580 continue
  1499. * le reste est nul
  1500. *
  1501. * remplissage du champoint de depi (jeu)
  1502. *
  1503. ipt7.num(1,nelch) = ipt1.num(1,iel1)
  1504. vpocha(nelch,1) = (xjeu - xjr) * SuPIn
  1505. vpocha(nelch,2) = SuPIn
  1506. IF (MELVA2.NE.0) THEN
  1507. NEL2 = min (iel1,NELA)
  1508. VPOCHA(nelch,3) = max(MELVA2.VELCHE(NPTELA,NEL2),0.D0)*SuPIn
  1509. ENDIF
  1510. C*DBG-F call prfaible(ipt1,iel,nPIn0,mfaible)
  1511. *
  1512. 520 CONTINUE
  1513. 519 CONTINUE
  1514. *
  1515. * Ajustement au plus juste puis desactivation des segments lies
  1516. * a la rigidite du type 0
  1517. if (nelri0.ne.nelrig) then
  1518. nelrig = nelri0
  1519. segadj,xmatri
  1520. nbelem = nelri0
  1521. nbnn = 7
  1522. segadj,meleme
  1523. endif
  1524. segdes,xmatri
  1525. *
  1526. ** write(ioimp,*) ' nb relation type 0 ',nelri0
  1527. * if (nelri0.eq.0) irigel(6,nrigel)=0
  1528. *
  1529. * Destruction des segments locaux devenus inutiles
  1530. segsup,mfaible
  1531. *
  1532. C* GOTO 1000
  1533.  
  1534. *=======================================================================
  1535. *= Fin du sous-programme :
  1536. *=======================================================================
  1537. 1000 CONTINUE
  1538. *
  1539. * Ajustement au plus juste du chpoint de depi (jeu) : mpoval et ipt7
  1540. * puis desactivation du chpoint
  1541. if (vpocha(/1).ne.nelch) then
  1542. n = nelch
  1543. nc=vpocha(/2)
  1544. segadj,mpoval
  1545. nbnn = 1
  1546. nbelem = nelch
  1547. nbsous = 0
  1548. nbref = 0
  1549. segadj,ipt7
  1550. ** write(6,*) ' ipt7 dans impo32 '
  1551. endif
  1552. * Desctivation de la matrice de raideur de contact
  1553. segdes,mrigid
  1554. *
  1555. * Reunion des relations portant sur le meme multiplicateur de lagrange
  1556. *
  1557. call impofu(MRIGID,MCHPOI)
  1558. * write(6,*) ' apres impofu dans impo32 '
  1559. * call ecchpo(mchpoi,1)
  1560. * call prrigi(mrigid,1)
  1561.  
  1562.  
  1563. *
  1564. * Nettoyage eventuel des termes petits dans les relations
  1565. *
  1566. * A voir plus tard. Frig3c suppose que la relation est complete.
  1567. *
  1568. ** call relasi(mrigid)
  1569. *
  1570. * Ecriture des resultats de l'operateur :
  1571. call ecrobj('RIGIDITE',mrigid)
  1572. call actobj('CHPOINT ',mchpoi,1)
  1573. call ecrobj('CHPOINT ',mchpoi)
  1574. *
  1575. 900 CONTINUE
  1576. end
  1577.  
  1578.  
  1579.  
  1580.  
  1581.  
  1582.  
  1583.  
  1584.  
  1585.  
  1586.  
  1587.  
  1588.  
  1589.  
  1590.  
  1591.  
  1592.  
  1593.  
  1594.  
  1595.  
  1596.  
  1597.  
  1598.  
  1599.  
  1600.  
  1601.  
  1602.  

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