Télécharger redcon.eso

Retour à la liste

Numérotation des lignes :

  1. C REDCON SOURCE CHAT 12/01/17 21:15:53 7248
  2. subroutine redcon(mmodel,mmode2)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. -INC CCOPTIO
  6. -INC SMCOORD
  7. -INC SMELEME
  8. -INC CCREEL
  9. -INC SMMODEL
  10. character*4 modepl(3),moforc(3)
  11. data modepl /'UX ','UY ','UZ '/
  12. data moforc /'FX ','FY ','FZ '/
  13. *
  14. * lecture du maillage support des conditions
  15. *
  16. segact mmodel
  17. segini,mmode2=mmodel
  18. imoin=0
  19. do 1 imo=1,kmodel(/1)
  20. imodel=kmodel(imo)
  21. segact imodel
  22. do ifo=1,formod(/2)
  23. if( formod(ifo).eq.'CONTACT' ) go to 2
  24. enddo
  25. mmode2.kmodel(imo-imoin)=imodel
  26. go to 1
  27. 2 continue
  28. * on va reduire le maillage
  29. ipt1=imamod
  30. segact ipt1
  31. segini,ipt2=ipt1
  32. igard=0
  33. nbelem=ipt1.num(/2)
  34. if(idim.eq.3) then
  35. *
  36. * boucle sur le maillage support en 3D
  37. *
  38.  
  39. do 10 iel=1,nbelem
  40. ip1=ipt1.num(2,iel)
  41. ip2=ipt1.num(3,iel)
  42. ip3=ipt1.num(4,iel)
  43. xp1=xcoor((ip1-1)*(idim+1)+1)
  44. yp1=xcoor((ip1-1)*(idim+1)+2)
  45. zp1=xcoor((ip1-1)*(idim+1)+3)
  46. xp2=xcoor((ip2-1)*(idim+1)+1)
  47. yp2=xcoor((ip2-1)*(idim+1)+2)
  48. zp2=xcoor((ip2-1)*(idim+1)+3)
  49. xp3=xcoor((ip3-1)*(idim+1)+1)
  50. yp3=xcoor((ip3-1)*(idim+1)+2)
  51. zp3=xcoor((ip3-1)*(idim+1)+3)
  52. *
  53. * verification si autre point dans la zone de selection
  54. jp=ipt1.num(5,iel)
  55. * write (6,*) 'mail ',ip1,ip2,ip3,jp
  56. xp =xcoor((jp-1)*(idim+1)+1)
  57. yp =xcoor((jp-1)*(idim+1)+2)
  58. zp =xcoor((jp-1)*(idim+1)+3)
  59. * verification que pas relation du point sur lui meme
  60. if (jp.eq.ip1) goto 10
  61. if (jp.eq.ip2) goto 10
  62. if (jp.eq.ip3) goto 10
  63. * verif distance au centre de gravite
  64. xg=(xp1+xp2+xp3)/3d0
  65. yg=(yp1+yp2+yp3)/3d0
  66. zg=(zp1+zp2+zp3)/3d0
  67. d1=sqrt((xg-xp1)**2+(yg-yp1)**2+(zg-zp1)**2)
  68. d2=sqrt((xg-xp2)**2+(yg-yp2)**2+(zg-zp2)**2)
  69. d3=sqrt((xg-xp3)**2+(yg-yp3)**2+(zg-zp3)**2)
  70. dp=sqrt((xg-xp)**2+(yg-yp)**2+(zg-zp)**2)/2d0
  71. if (dp.gt.d1.and.dp.gt.d2.and.dp.gt.d3) goto 10
  72. * write (6,*) 'contact test distance ok',dp,d1,d2,d3
  73. * verif position par rapport aux cotés
  74. * cote 1 2
  75. xv=xp2-xp1
  76. yv=yp2-yp1
  77. zv=zp2-zp1
  78. dv=sqrt((xv**2)+(yv**2)+(zv**2))
  79. if(.not.(dv.lt.1d-10).and.
  80. > .not.(dv.gt.1d-10))
  81. > write (6,*) 'prob 1 dans impo32'
  82. if (abs(dv).le.xpetit) write (6,*) 'prob 1 imp32'
  83. xv=xv/dv
  84. yv=yv/dv
  85. zv=zv/dv
  86. scal=(xp-xp1)*xv+(yp-yp1)*yv+(zp-zp1)*zv
  87. xm=xp1+scal*xv
  88. ym=yp1+scal*yv
  89. zm=zp1+scal*zv
  90. scal=(xg-xm)*xv+(yg-ym)*yv+(zg-zm)*zv
  91. xn=(xg-xm)-scal*xv
  92. yn=(yg-ym)-scal*yv
  93. zn=(zg-zm)-scal*zv
  94. dn=sqrt(xn**2+yn**2+zn**2)
  95. scal=(xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn
  96. dpm=sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2)
  97. if (dpm.gt.1d-3*dv.and.scal .lt.-0.2d0*dn*dpm) goto 10
  98. * write (6,*) 'contact test position 1 ok',scal/(dn*dpm),scal,
  99. * > dn,dpm
  100. * cote 2 3
  101. xv=xp3-xp2
  102. yv=yp3-yp2
  103. zv=zp3-zp2
  104. dv=sqrt((xv**2)+(yv**2)+(zv**2))
  105. if(.not.(dv.lt.1d-10).and.
  106. > .not.(dv.gt.1d-10))
  107. > write (6,*) 'prob 2 dans impo32'
  108. if (abs(dv).le.xpetit) write (6,*) 'prob 2 imp32'
  109. xv=xv/dv
  110. yv=yv/dv
  111. zv=zv/dv
  112. scal=(xp-xp2)*xv+(yp-yp2)*yv+(zp-zp2)*zv
  113. xm=xp2+scal*xv
  114. ym=yp2+scal*yv
  115. zm=zp2+scal*zv
  116. scal=(xg-xm)*xv+(yg-ym)*yv+(zg-zm)*zv
  117. xn=(xg-xm)-scal*xv
  118. yn=(yg-ym)-scal*yv
  119. zn=(zg-zm)-scal*zv
  120. dn=sqrt(xn**2+yn**2+zn**2)
  121. scal=(xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn
  122. dpm=sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2)
  123. if (dpm.gt.1d-3*dv.and.scal .lt.-0.2d0*dn*dpm) goto 10
  124. * write (6,*) 'contact test position 2 ok',scal/(dn*dpm),scal,
  125. * > dn,dpm
  126. * cote 3 1
  127. xv=xp1-xp3
  128. yv=yp1-yp3
  129. zv=zp1-zp3
  130. dv=sqrt((xv**2)+(yv**2)+(zv**2))
  131. if(.not.(dv.lt.1d-10).and.
  132. > .not.(dv.gt.1d-10))
  133. > write (6,*) 'prob 3 dans impo32'
  134. if (abs(dv).le.xpetit) write (6,*) 'prob 3 imp32'
  135. xv=xv/dv
  136. yv=yv/dv
  137. zv=zv/dv
  138. scal=(xp-xp3)*xv+(yp-yp3)*yv+(zp-zp3)*zv
  139. xm=xp3+scal*xv
  140. ym=yp3+scal*yv
  141. zm=zp3+scal*zv
  142. scal=(xg-xm)*xv+(yg-ym)*yv+(zg-zm)*zv
  143. xn=(xg-xm)-scal*xv
  144. yn=(yg-ym)-scal*yv
  145. zn=(zg-zm)-scal*zv
  146. dn=sqrt(xn**2+yn**2+zn**2)
  147. scal=(xp-xm)*xn+(yp-ym)*yn+(zp-zm)*zn
  148. dpm=sqrt((xp-xm)**2+(yp-ym)**2+(zp-zm)**2)
  149. if (dpm.gt.1d-3*dv.and.scal .lt.-0.2d0*dn*dpm) goto 10
  150. * write (6,*) 'contact test position 3 ok',scal/(dn*dpm),scal,
  151. * > dn,dpm
  152. * on a un point ou imposer la relation
  153. * Quelle est la relation ???
  154. * write (6,*) 'contact potentiel '
  155. *
  156. * direction de la relation
  157. *
  158.  
  159. xn=(yp1-yp2)*(zp2-zp3)-(zp1-zp2)*(yp2-yp3)
  160. yn=(zp1-zp2)*(xp2-xp3)-(xp1-xp2)*(zp2-zp3)
  161. zn=(xp1-xp2)*(yp2-yp3)-(yp1-yp2)*(xp2-xp3)
  162. dn=sqrt(xn**2+yn**2+zn**2)
  163. if(.not.(dn.lt.1d-10).and.
  164. > .not.(dn.gt.1d-10))
  165. > write (6,*) 'prob 4 dans impo32'
  166. if (abs(dn).le.xpetit) write (6,*) 'prob 4 imp32'
  167. xn=xn/dn
  168. yn=yn/dn
  169. zn=zn/dn
  170. if (abs(xn).lt.1d-10) xn=0.d0
  171. if (abs(yn).lt.1d-10) yn=0.d0
  172. if (abs(zn).lt.1d-10) zn=0.d0
  173. *
  174. * recherche du point a mettre en relation (coordonnees barycentriques)
  175. *
  176. scal=(xp-xg)*xn+(yp-yg)*yn+(zp-zg)*zn
  177. *
  178. * verif bon cote (a la tolerance pres)
  179. *
  180. if (scal.lt.-0.2d0*d1.and.scal.lt.-0.2d0*d2.and.
  181. > scal.lt.-0.2d0*d3) goto 10
  182.  
  183. *
  184. * on range dans le meleme
  185. *
  186.  
  187. igard=igard+1
  188. ipt2.num(1,igard)=ipt1.num(1,iel)
  189. ipt2.num(2,igard)=ip1
  190. ipt2.num(3,igard)=ip2
  191. ipt2.num(4,igard)=ip3
  192. ipt2.num(5,igard)=jp
  193. 10 continue
  194. segdes imodel,ipt1
  195. if(igard.eq.nbelem) then
  196. mmode2.kmodel(imo-imoin)=imodel
  197. segsup ipt2
  198. go to 1
  199. endif
  200. if( igard.eq.0) then
  201. imoin=imoin+1
  202. segsup ipt2
  203. go to 1
  204. endif
  205. nbelem=igard
  206. nbnn=ipt2.num(/1)
  207. nbref=0
  208. nbsous=0
  209. segadj ipt2
  210. segini,imode2=imodel
  211. imode2.imamod=ipt2
  212. kmodel(imo-imoin)=imode2
  213. segdes imode2,ipt2
  214. go to 1
  215. elseif(idim.eq.2) then
  216. if( ipt1.num(/1).eq.5) then
  217. * formulation faible
  218. do 210 iel=1,ipt1.num(/2)
  219. ip1=ipt1.num(2,iel)
  220. ip2=ipt1.num(3,iel)
  221. ip3=ipt1.num(4,iel)
  222. ip4=ipt1.num(5,iel)
  223. xp1=xcoor((ip1-1)*(idim+1)+1)
  224. yp1=xcoor((ip1-1)*(idim+1)+2)
  225. xp2=xcoor((ip2-1)*(idim+1)+1)
  226. yp2=xcoor((ip2-1)*(idim+1)+2)
  227. xp3=xcoor((ip3-1)*(idim+1)+1)
  228. yp3=xcoor((ip3-1)*(idim+1)+2)
  229. xp4=xcoor((ip4-1)*(idim+1)+1)
  230. yp4=xcoor((ip4-1)*(idim+1)+2)
  231. *
  232. * critere d'acceptation de l'éléments:
  233. * une extrémité d'un segment dans le cercle de sélection de l'autre segment
  234. * orientations respectives correctes des 2 segments
  235. *
  236. * verification (provisoire) que pas de noeuds doubles.
  237. if (ip1.eq.ip3) goto 210
  238. if (ip1.eq.ip4) goto 210
  239. if (ip2.eq.ip3) goto 210
  240. if (ip2.eq.ip4) goto 210
  241. *
  242. xl13=sqrt((xp3-xp1)**2+(yp3-yp1)**2)
  243. xl14=sqrt((xp4-xp1)**2+(yp4-yp1)**2)
  244. xl23=sqrt((xp3-xp2)**2+(yp3-yp2)**2)
  245. xl24=sqrt((xp4-xp2)**2+(yp4-yp2)**2)
  246. sca312=(xp3-xp1)*(xp3-xp2)+(yp3-yp1)*(yp3-yp2)
  247. sca412=(xp4-xp1)*(xp4-xp2)+(yp4-yp1)*(yp4-yp2)
  248. sca134=(xp1-xp3)*(xp1-xp4)+(yp1-yp3)*(yp1-yp4)
  249. sca234=(xp2-xp3)*(xp2-xp4)+(yp2-yp3)*(yp2-yp4)
  250. if (sca312/(xl13*xl23).gt.0.7.and.
  251. > sca412/(xl14*xl24).gt.0.7.and.
  252. > sca134/(xl13*xl14).gt.0.7.and.
  253. > sca234/(xl23*xl24).gt.0.7) goto 210
  254. *
  255. scal=(xp1-xp2)*(xp3-xp4)+(yp1-yp2)*(yp3-yp4)
  256. if (scal.gt.0.) goto 210
  257.  
  258. * on a un element ou imposer la relation
  259. * Quelle est la relation ???
  260. *
  261. * direction du contact
  262. *
  263. xr=(xp2-xp1+xp3-xp4)
  264. yr=(yp2-yp1+yp3-yp4)
  265. sr2 = xr**2 + yr**2
  266. sr = sqrt (sr2)
  267. xr=xr/sr
  268. yr=yr/sr
  269. *
  270. * projection des points sur la direction du contact
  271. *
  272. xl1=xp1*xr+yp1*yr
  273. xl2=xp2*xr+yp2*yr
  274. xl3=xp3*xr+yp3*yr
  275. xl4=xp4*xr+yp4*yr
  276. *
  277. * distance des points à leurs projections
  278. *
  279. d1=(xp1-xl1*xr)*yr-(yp1-xl1*yr)*xr
  280. d2=(xp2-xl2*xr)*yr-(yp2-xl2*yr)*xr
  281. d3=(xp3-xl3*xr)*yr-(yp3-xl3*yr)*xr
  282. d4=(xp4-xl4*xr)*yr-(yp4-xl4*yr)*xr
  283. *
  284. * calcul de l'intersection des projections.
  285. *
  286. xm1=min(xl1,xl2)
  287. xm2=max(xl1,xl2)
  288. xm3=min(xl3,xl4)
  289. xm4=max(xl3,xl4)
  290. *
  291. * critere
  292. *
  293. xcr=max(xm2-xm1,xm4-xm3)*1d-4
  294. xi=max(xm1,xm3)
  295. xj=min(xm2,xm4)
  296. *
  297. if (xi.ge.xj-xcr) goto 210
  298.  
  299. * on range dans le meleme
  300. *
  301. igard=igard+1
  302. ipt2.num(1,igard)=ipt1.num(1,iel)
  303. ipt2.num(2,igard)=ip1
  304. ipt2.num(3,igard)=ip2
  305. ipt2.num(4,igard)=ip3
  306. ipt2.num(5,igard)=ip4
  307. *
  308. 210 continue
  309. segdes imodel,ipt1
  310. if(igard.eq.nbelem) then
  311. mmode2.kmodel(imo-imoin)=imodel
  312. segsup ipt2
  313. go to 1
  314. endif
  315. if( igard.eq.0) then
  316. imoin=imoin+1
  317. segsup ipt2
  318. go to 1
  319. endif
  320. nbelem=igard
  321. nbnn=ipt2.num(/1)
  322. nbref=0
  323. nbsous=0
  324. segadj ipt2
  325. segini,imode2=imodel
  326. imode2.imamod=ipt2
  327. kmodel(imo-imoin)=imode2
  328. segdes imode2,ipt2
  329. go to 1
  330. elseif(ipt1.num(/1).eq.4) then
  331. * formulation 1 point - 1 segment
  332. do 110 iel=1,ipt1.num(/2)
  333. ip1=ipt1.num(2,iel)
  334. ip2=ipt1.num(3,iel)
  335. xp1=xcoor((ip1-1)*(idim+1)+1)
  336. yp1=xcoor((ip1-1)*(idim+1)+2)
  337. xp2=xcoor((ip2-1)*(idim+1)+1)
  338. yp2=xcoor((ip2-1)*(idim+1)+2)
  339. *
  340. * verification si autre point dans le cercle de selection
  341. jp=ipt1.num(4,iel)
  342. * verification que pas relation du point sur lui meme
  343. if (jp.eq.ip1) goto 110
  344. if (jp.eq.ip2) goto 110
  345. xp=xcoor((jp-1)*(idim+1)+1)
  346. yp=xcoor((jp-1)*(idim+1)+2)
  347. scal=(xp-xp1)*(xp-xp2)+(yp-yp1)*(yp-yp2)
  348. xl1=sqrt((xp-xp1)**2+(yp-yp1)**2)
  349. xl2=sqrt((xp-xp2)**2+(yp-yp2)**2)
  350. scal=scal/(xl1*xl2)
  351. if (scal.gt.0.7) goto 110
  352. igard=igard+1
  353. ipt2.num(1,igard)=ipt1.num(1,iel)
  354. ipt2.num(2,igard)=ip1
  355. ipt2.num(3,igard)=ip2
  356. ipt2.num(4,igard)=jp
  357. 110 continue
  358. segdes imodel,ipt1
  359. if(igard.eq.nbelem) then
  360. mmode2.kmodel(imo-imoin)=imodel
  361. segsup ipt2
  362. go to 1
  363. endif
  364. if( igard.eq.0) then
  365. imoin=imoin+1
  366. segsup ipt2
  367. go to 1
  368. endif
  369. nbelem=igard
  370. nbnn=ipt2.num(/1)
  371. nbref=0
  372. nbsous=0
  373. segadj ipt2
  374. segini,imode2=imodel
  375. imode2.imamod=ipt2
  376. mmode2.kmodel(imo-imoin)=imode2
  377. segdes imode2,ipt2
  378. go to 1
  379. endif
  380. endif
  381. 1 continue
  382. segdes mmodel
  383. if(imoin.ne.0) then
  384. n1=kmodel(/1)-imoin
  385. if(n1.eq.0) then
  386. call erreur (21)
  387. return
  388. endif
  389. segadj mmode2
  390. endif
  391. segdes mmode2
  392. return
  393. end
  394.  
  395.  

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