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

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