Télécharger reltuy.eso

Retour à la liste

Numérotation des lignes :

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

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