Télécharger reltuy.eso

Retour à la liste

Numérotation des lignes :

  1. C RELTUY SOURCE CHAT 12/09/11 21:15:11 7499
  2. subroutine reltuy
  3. implicit real*8(a-h,o-z)
  4. implicit integer (i-n)
  5. -INC CCOPTIO
  6. -INC SMRIGID
  7. -INC SMELEME
  8. -INC SMCOORD
  9. -INC CCREEL
  10. segment icpr(xcoor(/1)/(idim+1))
  11. SEGMENT ISEG3
  12. INTEGER NIZO(NxyZo+1)
  13. ENDSEGMENT
  14. * SEGMENT ISEG4
  15. * INTEGER NUMZO(NxyZo)
  16. * ENDSEGMENT
  17. SEGMENT ISEG5
  18. INTEGER NNMEL(ILON),IDEJ(nel),ipos(na+1)
  19. ENDSEGMENT
  20. SEGMENT ISEG6
  21. INTEGER numpo(illo)
  22. ENDSEGMENT
  23. segment izon(na)
  24. DATA EPSI/1.D-4/
  25.  
  26. *
  27. * lecture des donnees on met le maillage de seg2 ou seg3 dans ipt1
  28. *
  29. call lirobj( 'MAILLAGE',ipt1,1,iretou)
  30. if(ierr.ne.0) return
  31. call lirobj('MAILLAGE' ,ipt2,1,iretou)
  32. if(ierr.ne.0) return
  33. segact ipt1,ipt2
  34. if(ipt1.lisous(/1).ne.0) then
  35. if(ipt2.itypel.ne.2.and.ipt2.itypel.ne.3)then
  36. call erreur (19)
  37. return
  38. endif
  39. ipta=ipt2
  40. ipt2=ipt1
  41. ipt1=ipta
  42. else
  43. if (ipt1.itypel.ne.2.and.ipt1.itypel.ne.3) then
  44. if(ipt2.itypel.ne.2.and.ipt2.itypel.ne.3)then
  45. call erreur(19)
  46. return
  47. else
  48. ipta=ipt2
  49. ipt2=ipt1
  50. ipt1=ipta
  51. endif
  52. endif
  53. endif
  54. call change (ipt2,1)
  55. if( ipt1.itypel.eq.3) call change (ipt1,2)
  56. if( ipt1.lisous(/1).ne.0) then
  57. * on suppose un seul type d'elements pour la tuyauterir
  58. call erreur(19)
  59. return
  60. endif
  61. npipt2=ipt2.num(/2)
  62. *
  63. *
  64. * lecture facultative de la plus grande section d'un tuyau
  65. *
  66. call lirree(sect,1,iretou)
  67. if(ierr.ne.0) return
  68. *
  69. *pour zoner on prend la longueurn max d'un element de ajouter a sect
  70. *
  71. xlon=0.d0
  72. ider=ipt1.num(/1)
  73. do iel=1,ipt1.num(/2)
  74. ia= ipt1.num(1,iel)
  75. ib=ipt1.num(ider,iel)
  76. ia1=ia*(idim+1)-idim
  77. xa=xcoor(ia1)
  78. ya=xcoor ( ia1+1)
  79. za=xcoor ( ia1+2)
  80. ib1=ib*(idim+1)-idim
  81. xb=xcoor(ib1)
  82. yb=xcoor ( ib1+1)
  83. zb=xcoor ( ib1+2)
  84. xlo= sqrt( (xa-xb)*(xa-xb)+(ya-yb)*(ya-yb)+(za-zb)*(za-zb))
  85. if( xlo . gt. xlon) xlon=xlo
  86. enddo
  87. xpre=epsi*xlon
  88. xcri= sqrt(xlon*xlon + sect*sect)
  89. *
  90. *
  91. * zonage on cherche les min et max suivant les trois directions et la taille
  92. *max des elements de la ligne de tuyauterie suivant les trois directions
  93. *
  94. *
  95. ITR = 0
  96. XZO=0.D0
  97. YZO=0.D0
  98. ZZO=0.D0
  99. XZA=XGRAND
  100. YZA=XGRAND
  101. ZZA=XGRAND
  102. XMI=XGRAND
  103. XMA=-XGRAND
  104. YMI=XGRAND
  105. YMA=-XGRAND
  106. if(idim.eq.3) then
  107. ZMI=XGRAND
  108. ZMA=-XGRAND
  109. else
  110. zmi=0.d0
  111. zma=0.d0
  112. endif
  113. meleme=ipt1
  114. ipasag=0
  115. idim1=idim+1
  116. 254 ipasag=ipasag+1
  117. if(ipasag.eq.1) then
  118. meleme=ipt1
  119. else
  120. meleme=ipt2
  121. endif
  122. nbnn= num(/1)
  123. nbel=num(/2)
  124. do 200 i1=1,nbel
  125. DO 201 I2 = 1,NBNN
  126. IB=NUM(I2,I1)
  127. IA=(IB-1)*IDIM1
  128. IF(XCOOR(IA+1).LT.XMI) XMI=XCOOR(IA+1)
  129. IF(XCOOR(IA+1).GT.XMA) XMA=XCOOR(IA+1)
  130. IF(XCOOR(IA+2).LT.YMI) YMI=XCOOR(IA+2)
  131. IF(XCOOR(IA+2).GT.YMA) YMA=XCOOR(IA+2)
  132. IF ( IDIM.EQ.3 ) THEN
  133. IF(XCOOR(IA+3).LT.ZMI) ZMI=XCOOR(IA+3)
  134. IF(XCOOR(IA+3).GT.ZMA) ZMA=XCOOR(IA+3)
  135. ENDIF
  136. 201 CONTINUE
  137. 200 CONTINUE
  138. If(ipasag.eq.1) go to 254
  139. nzo=1
  140. xmi= xmi-( xma-xmi ) / 1.d6
  141. ymi= ymi-( yma-ymi ) / 1.d6
  142. xma=xma + ( xma-xmi ) / 1.d6
  143. yma=yma + ( yma-ymi ) / 1.d6
  144. nxo= int(( xma - xmi ) / xcri) +1
  145. nyo= int(( yma - ymi ) / xcri) + 1
  146. if(idim.eq.3) then
  147. zmi= zmi-( zma-zmi ) / 1.d6
  148. zma=zma+( zma-zmi ) / 1.d6
  149. nzo= int (( zma - zmi ) / xcri) + 1
  150. endif
  151. nxyo=nxo*nyo
  152. nxyzo= nxo*nyo*nzo
  153. * write(6,*) ' xcri nxo nyo nzo nxyo nxyzo'
  154. * write(6,*) xcri, nxo ,nyo, nzo, nxyo ,nxyzo
  155. * write(6,*) ' xmi xma ymi yma zmi zma',xmi,xma,ymi,yma,zmi,zma
  156. *
  157. * on fait une numerotation locale pour les points de la ligne
  158. * on ne s'occupe que des points extremes
  159. *
  160. na=0
  161. segini icpr
  162. do iel=1,ipt1.num(/2)
  163. do ip=1,ipt1.num(/1)
  164. ia= ipt1.num(ip,iel)
  165. if(icpr(ia).eq.0) then
  166. na=na+1
  167. icpr(ia)=na
  168. endif
  169. enddo
  170. enddo
  171. * write(6,*) ' na nxyo nxo' , na,nxyo , nxo
  172. *
  173. * on cherche dans quelle zone est chacun des points de la tuyauterie
  174. * et on dresse le tableau la liste des points par zone
  175.  
  176. segini iseg3,izon
  177. zp=0.d0
  178. do iu=1,icpr(/1)
  179. if(icpr(iu).ne.0) then
  180. ip=icpr(iu)
  181. xp= xcoor( (iu-1)*(idim+1) +1)
  182. yp=xcoor( (iu-1)*(idim+1) +2)
  183. if( idim.eq.3)ZP= xcoor((iu-1)*(idim+1) +3)
  184. nn= (int(( zp-zmi)/xcri) )*nxyo
  185. $ +(int(( yp-ymi)/xcri) )*nxo
  186. $ +(int(( xp-xmi)/xcri) + 1)
  187. nizo(nn)=nizo(nn)+1
  188. izon(ip)=nn
  189. * write(6,*) 'xp yp zp iu ip izon(ip) ',xp,yp,zp ,iu,ip,izon(ip)
  190. endif
  191. enddo
  192. do ia=2,nxyzo+1
  193. nizo(ia)=nizo(ia)+nizo(ia-1)
  194. enddo
  195. * write(6,*) ' nizo ' , (nizo(iou),iou=1,nizo(/1))
  196. illo=nizo(nxyzo+1)
  197. segini iseg6
  198. do ip1=1,icpr(/1)
  199. inp=icpr(ip1)
  200. if( inp.ne.0) then
  201. nn=izon(inp)
  202. ipla=nizo(nn)
  203. nizo(nn)=nizo(nn)-1
  204. numpo(ipla)=ip1
  205. endif
  206. enddo
  207. * write(6,*) ' nizo ' , (nizo(iou),iou=1,nizo(/1))
  208. * write(6,*) ' numpo ',(numpo(iou),iou=1,numpo(/1))
  209. * on construit le tableua donnant le numero des elements touchant un noeud
  210. * en numerotation locale
  211. ilon=2*na
  212. nel=ipt1.num(/2)
  213. segini iseg5
  214. do ia=1,ipt1.num(/2)
  215. do ib=1,ipt1.num(/1)
  216. ip1=ipt1.num(ib,ia)
  217. iy=icpr(ip1)
  218. ipos(iy)=ipos(iy)+1
  219. enddo
  220. enddo
  221. do ia=2,na+1
  222. ipos(ia)=ipos(ia-1)+ipos(ia)
  223. enddo
  224. do ia=1,ipt1.num(/2)
  225. do ib=1,ipt1.num(/1)
  226. ip1=ipt1.num(ib,ia)
  227. iy=icpr(ip1)
  228. ipla=ipos(iy)
  229. ipos(iy)=ipos(iy)-1
  230. nnmel(ipla)=ia
  231. enddo
  232. enddo
  233. * write(6,*) ' ipos ' , ( ipos(iou),iou=1,ipos(/1))
  234. * write(6,*) ' nnmel ' , ( nnmel(iou),iou=1,nnmel(/1))
  235. *
  236. * on cree tout de suite les points des multiplicateur
  237. *
  238. iposmu=xcoor(/1)/(idim+1)
  239. nbpts=iposmu + npipt2
  240. segadj mcoord
  241.  
  242. *
  243. * pour les relations il peut y avoir des relations sur tout le segment ou sur
  244. * seulement un point. ON cree les deux types, on ajustera plus tard.
  245. *
  246. nbnn= ipt1.num(/1)+2
  247. nbelem= ipt2.num(/2)
  248. nbsous=0
  249. nbref=0
  250. segini ipt4
  251. ipt4.itypel=22
  252. nrigel=2
  253. segini mrigid
  254. mtymat='RIGIDITE'
  255. coerig(1)=1.d0
  256. coerig(2)=1.d0
  257. irigel(1,1)=ipt4
  258. nligrp=nbnn
  259. nligrd=nbnn
  260. segini descr
  261. irigel(3,1)=descr
  262. lisinc(1)='LX'
  263. lisdua(1)='FLX'
  264. noelep(1)=1
  265. noeleD(1)=1
  266. do iu=2,nbnn
  267. noelep(iu)=iu
  268. noeled(iu)=iu
  269. lisinc(iu)='T'
  270. lisdua(iu)='Q'
  271. enddo
  272. segdes descr
  273. nelrig= nbelem
  274. segini xmatr4
  275. irigel(4,1)= xmatr4
  276. irigel(5,1)=nifour
  277. nel4=0
  278. nel5=0
  279. nbnn=3
  280. * write(6,*) ' à la creation de ipt5 nbelem ' , nbelem
  281. segini ipt5
  282. ipt5.itypel=22
  283. irigel(1,2)=ipt5
  284. nligrp=nbnn
  285. nligrd=nbnn
  286. segini xmatr5
  287. irigel(4,2)=xmatr5
  288. irigel(5,2)=nifour
  289. segini descr
  290. irigel(3,2)=descr
  291. lisinc(1)='LX'
  292. lisdua(1)='FLX'
  293. noelep(1)=1
  294. noeleD(1)=1
  295. noelep(2)=2
  296. noeled(2)=2
  297. noelep(3)=3
  298. noeled(3)=3
  299. lisinc(2)='T'
  300. lisinc(3)='T'
  301. lisdua(2)='Q'
  302. lisdua(3)='Q'
  303. segdes descr
  304. nbnn=ipt1.num(/1)
  305. nz=1
  306. do ipp=1,ipt2.num(/2)
  307. ip= ipt2.num(1,ipp)
  308. iqv=0
  309. xp= xcoor( (ip-1)*(idim+1) +1)
  310. yp=xcoor( (ip-1)*(idim+1) +2)
  311. zp=0.d0
  312. if(idim.eq.3)ZP= xcoor((ip-1)*(idim+1) +3)
  313. nx= int(( xp-xmi)/xcri) + 1
  314. ny= int(( yp-ymi)/xcri) + 1
  315. nz= int(( zp-zmi)/xcri) + 1
  316. nn= (nz-1)*nxyo +(ny-1)*nxo + nx
  317. * write(6,*) ' point ' , ip , ' zonze ' , nn, xp, yp,zp,nzo
  318. * write(6,*) ' point , ip ,nx ny nz ' , nx , ny , nz
  319. ield=0
  320. xdi=xgrand
  321. nnxmi= max ( nx-1,1)
  322. nnxma= min (nx+1,nxo)
  323. nnymi= max(1,ny-1)
  324. nnyma= min(nyo,nY+1)
  325. nnzmi= max(1,nz-1)
  326. nnzma= min(nzo,nz+1)
  327. * write(6,*) 'mima', nnxmi,nnxma,nnymi,nnyma,nnzmi,nnzma
  328. do ix=nnxmi,nnxma
  329. do iy=nnymi,nnyma
  330. do iz=nnzmi,nnzma
  331. nnes=(iz-1)*nxyo + (iy-1)*nxo+ix
  332. itp =nizo(nnes+1)-nizo(nnes)
  333. * write(6,* ) ' ip ix iy iz nnes itp' , ix,iy,iz, nnes,itp
  334. if (itp.ne.0) then
  335. idec= nizo(nnes)
  336. do im=1,itp
  337. ippp=numpo( idec + im)
  338. ipl=icpr(ippp)
  339. * write(6,*) ' ippp ipl ' , ippp,ipl
  340. iplael=ipos(ipl)+1
  341. do ielt=iplael,ipos(ipl+1)
  342. iel= nnmel(ielt)
  343. * write(6,*) ' ielt iel nbnn' , ielt,iel,nbnn
  344. ipa=(ipt1.num(1,iel)-1)*(idim+1)
  345. ipb=(ipt1.num(nbnn,iel)-1)*(idim+1)
  346. * write(6,*) ' ipa ipb ' , ipa/(idim+1), ipb/(idim+1)
  347. xa = xcoor(ipa+1)
  348. ya = xcoor(ipa+2)
  349. xb = xcoor(ipb+1)
  350. yb = xcoor(ipb+2)
  351. if(idim.eq.3) then
  352. za = xcoor(ipa+3)
  353. zb = xcoor(ipb+3)
  354. else
  355. za=0.d0
  356. zb=0.d0
  357. endif
  358. *
  359. * calcul de la distance de ip a la droit A B et verif si pt interne ou non
  360. *
  361. xab=xb-xa
  362. yab=yb-ya
  363. zab=zb-za
  364. xlo= sqrt(xab*xab + yab * yab + zab * zab)
  365. xab=xab/xlo
  366. yab=yab/xlo
  367. zab=zab/xlo
  368. xloi= xab*( xp-xa)+yab*(yp-ya)+zab*(zp-za)
  369. xloi=xloi/xlo
  370. * if( ip.eq.118476) then
  371. * write(6,*) ' xloi xpre xlo ' , xloi, xpre,xlo
  372. * write(6,*)' ipa ipb ',ipt1.num(1,iel),ipt1.num(nbnn,iel)
  373. * write(6,*) ' xab yab zab ',xab,yab,zab
  374. * write(6,*) ' xa ya za xp yp zp',xa,ya,za,xp,yp,zp
  375. * endif
  376. if( xloi.ge.epsi.and. xloi.le.(1.D0-epsi) ) then
  377. vx= yab*(zp-za) - zab*(yp-ya)
  378. vy= zab*(xp-xa) - xab*(zp-za)
  379. vz= xab*(yp-ya) - yab*(xp-xa)
  380. xdis=sqrt(vx*vx + vy*vy + vz*vz)
  381. * if( ip.eq.118476) then
  382. * write(6,*) ' on a trouve un candidat xdis' , xdis
  383. * endif
  384. if( xdis.lt.xdi) then
  385. ield=iel
  386. xdi=xdis
  387. xrap=xloi
  388. endif
  389. endif
  390. enddo
  391. enddo
  392. endif
  393. enddo
  394. enddo
  395. enddo
  396. * a-t-on trouvé un candidat?
  397. if(ield.ne.0) then
  398. iposmu=iposmu+1
  399. nel4=nel4+1
  400. ipt4.num(1,nel4)=iposmu
  401. ipt4.num(2,nel4)=ip
  402. ipt4.num(3,nel4)=ipt1.num(1,ield)
  403. ipt4.num(4,nel4)=ipt1.num(2,ield)
  404. xmatr4.re(1,2,nel4)=1.d0
  405. xmatr4.re(1,3,nel4)= xrap -1.D0
  406. xmatr4.re(1,4,nel4)= -xrap
  407. xmatr4.re(2,1,nel4)=1.d0
  408. xmatr4.re(3,1,nel4)= xrap -1.D0
  409. xmatr4.re(4,1,nel4)= -xrap
  410. else
  411. * il faut recommencer pour trouver le point le plus proche et ecrire une
  412. * relation sur un seul point
  413. xdi=xgrand
  414. * write(6,*) 'nnxmi,nnxma',nnxmi,nnxma
  415. * write(6,*) 'nnymi,nnyma',nnymi,nnyma
  416. * write(6,*) 'nnzmi,nnzma',nnzmi,nnzma
  417. do ix=nnxmi,nnxma
  418. do iy=nnymi,nnyma
  419. do iz=nnzmi,nnzma
  420. nnes= (iz-1)*nxyo + (iy-1)*nxo+ix
  421. * write(6,*) 'nizo(nnes)+1,nizo(nnes+1)',nizo(nnes),nizo(nnes+1)
  422. do iqz=nizo(nnes)+1,nizo(nnes+1)
  423. iq=numpo(iqz)
  424.  
  425. iqq=(iq-1)*(idim+1)
  426. xq=xcoor(iqq+1)
  427. yq=xcoor(iqq+2)
  428. zq=0.d0
  429. if( idim.eq.3) zq=xcoor(iqq+3)
  430. xdis=(xp-xq)*(xp-xq) + (yp-yq)*(yp-yq) + ( zp-zq)*(zp-zq)
  431. * write(6,*) ' point essaye iq xdis ' , iq ,xdis , xdi
  432. if( xdi.gt.xdis) then
  433. xdi=xdis
  434. iqv=iq
  435. endif
  436. enddo
  437. enddo
  438. enddo
  439. enddo
  440. nel5=nel5+1
  441. iposmu=iposmu+1
  442. ipt5.num(1,nel5)=iposmu
  443. ipt5.num(2,nel5)=ip
  444. ipt5.num(3,nel5)=iqv
  445. xmatr5.re(1,2,nel5)=1.d0
  446. xmatr5.re(2,1,nel5)=1.D0
  447. xmatr5.re(1,3,nel5)=-1.d0
  448. xmatr5.re(3,1,nel5)=-1.d0
  449. endif
  450. enddo
  451. *
  452. * ajustement des objets rigidité
  453. *
  454.  
  455. * write(6,*) ' nbelem nel4 nel5 ' , nbelem,nel4 , nel5
  456. * write(6,*) ' irigel', ( irigel(iou,1),iou=1,8)
  457. if(nel4.ne.0) then
  458. if(ipt4.num(/2).ne.nel4) then
  459. nelrig=nel4
  460. nbelem=nel4
  461. nbnn= ipt4.num(/1)
  462. segadjipt4
  463. nligrp=xmatr4.re(/1)
  464. nligrd=xmatr4.re(/2)
  465. segadj xmatr4
  466. nelrig= nel5
  467. nbelem=nel5
  468. nbnn=ipt5.num(/1)
  469. segadj ipt5
  470. nligrp=xmatr5.re(/1)
  471. nligrd=xmatr5.re(/2)
  472. segadj xmatr5
  473. else
  474. * write(6,*) ' on passe bien la'
  475. nrigel=1
  476. segadj mrigid
  477. * write(6,*) ' irigel', ( irigel(iou,1),iou=1,8)
  478. endif
  479. else
  480. do io=1,irigel(/1)
  481. irigel(io,1)=irigel(io,2)
  482. enddo
  483. nrigel=1
  484. segadj mrigid
  485. endif
  486. * write(6,*) ' irigel', ( irigel(iou,1),iou=1,8)
  487. segdes ipt4,ipt5,xmatr4,xmatr5,mrigid
  488. segsup izon,iseg3,iseg5,iseg6,icpr
  489. call ecrobj('RIGIDITE',mrigid)
  490. return
  491. end
  492.  
  493.  
  494.  

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