Télécharger proob4.eso

Retour à la liste

Numérotation des lignes :

  1. C PROOB4 SOURCE CHAT 05/01/13 02:34:16 5004
  2. C CE PROGRAMME PERMET LA PROJECTION CONIQUE D'UN MAILLAGE 3D SUR
  3. C UNE SUFACE CARACTERISEE PAR UN AUTRE MAILLAGE 3D.
  4.  
  5. C ENTREES/ MAI2 : UN MAILLAGE DE TYPE MELEME A PROJETER.
  6. C MAI1 : UN MAILLAGE REPRESENTANT LA SURFACE SUR LAQUELLE MAI2
  7. C EST PROJETE.
  8. C IP1 : UN POINT DEFINISSANT L'OEIL DE LA PROJECTION.
  9.  
  10. C DATE: SEPTEMBRE 1996.
  11.  
  12. C **********************************************************************
  13. Subroutine proob4(mai2,mai1,ip1)
  14. IMPLICIT INTEGER(I-N)
  15. IMPLICIT REAL*8(A-H,O-Z)
  16.  
  17. -INC CCOPTIO
  18. -INC SMELEME
  19. -INC SMCOORD
  20. -INC CCREEL
  21. *-
  22. SEGMENT ISEG1
  23. REAL*8 XLIM(2,NBEL),YLIM(2,NBEL),ZLIM(2,NBEL)
  24. ENDSEGMENT
  25.  
  26. SEGMENT ISEG3
  27. INTEGER NIZO(NZO+1)
  28. ENDSEGMENT
  29.  
  30. SEGMENT ISEG4
  31. INTEGER NUMZO(NZO)
  32. ENDSEGMENT
  33.  
  34. SEGMENT ISEG5
  35. INTEGER NNMEL(ILON),IDEJ(NZO)
  36. ENDSEGMENT
  37.  
  38. SEGMENT ISEG6
  39. REAL*8 AM(4,4),AL(4),A(4),B(4),C(4)
  40. REAL*8 XPU(2),P(3)
  41. ENDSEGMENT
  42.  
  43. SEGMENT MCOOR1
  44. REAL*8 XCOOR1(XCOOR(/1))
  45. ENDSEGMENT
  46. C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE.
  47.  
  48. SEGMENT MCOOR2
  49. REAL*8 XCOOR2(XCOOR(/1))
  50. ENDSEGMENT
  51. C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE.
  52.  
  53. SEGMENT MCOOR3
  54. REAL*8 XCOOR3(XCOOR(/1))
  55. ENDSEGMENT
  56. C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE.
  57.  
  58. SEGMENT MCOOR4
  59. REAL*8 XCOOR4(XCOOR(/1))
  60. ENDSEGMENT
  61. C LA DERNIERE COORDONNEE D'UN POINT N'EST PAS UTILISEE.
  62.  
  63. SEGMENT MCOR
  64. INTEGER ICOR(IONMAX+1)
  65. ENDSEGMENT
  66.  
  67. SEGMENT MCORRE
  68. INTEGER ICORRE((XCOOR(/1)/(IDIM+1)))
  69. ENDSEGMENT
  70.  
  71. SEGMENT ICP1
  72. INTEGER ICPR1(XCOOR(/1)/(IDIM+1))
  73. ENDSEGMENT
  74.  
  75. SEGMENT ICP2
  76. INTEGER ICPR2(XCOOR(/1)/(IDIM+1))
  77. ENDSEGMENT
  78.  
  79. SEGMENT MATRIC
  80. REAL*8 TABMAT(3,3),V(3),W(3),Q(3),invma(3,3)
  81. ENDSEGMENT
  82.  
  83. SEGMENT CRIT
  84. INTEGER ELECRI(NBELEM),TABCRIT(NBELEM+1)
  85. ENDSEGMENT
  86.  
  87.  
  88. segini matric
  89. nbpts=xcoor(/1)/(idim+1)
  90. zchoix=100.
  91. w(1)=xcoor((ip1-1)*(idim+1)+1)
  92. w(2)=xcoor((ip1-1)*(idim+1)+2)
  93. w(3)=xcoor((ip1-1)*(idim+1)+3)
  94.  
  95. C PASSAGE DU MAILLAGE MAI1 EN TRI3
  96. call change(mai1,4)
  97.  
  98. C LISTAGE DES POINTS APPARTENANT A MAI1
  99. segini icp1
  100. nbpts = xcoor(/1)/( idim+1)
  101. meleme=mai1
  102. segact meleme
  103. nbelem=num(/2)
  104. nbnn = num(/1)
  105. Do 10 i=1,nbelem
  106. Do 20 j=1,nbnn
  107. ia=num(j,i)
  108. icpr1(ia)=1
  109. 20 continue
  110. 10 continue
  111.  
  112. C LISTAGE DES POINTS APPARTENANT A MAI2
  113. ipt1=mai2
  114. segact ipt1
  115. segini icp2
  116. nbelem= ipt1.num(/2)
  117. nbnn = ipt1.num(/1)
  118. Do 30 i=1,nbelem
  119. Do 40 j=1,nbnn
  120. ia=ipt1.num(j,i)
  121. If (icpr1(ia).eq.1) then
  122. icpr2(ia)=2
  123. else
  124. icpr2(ia)=1
  125. Endif
  126. 40 continue
  127. 30 continue
  128.  
  129. C POSITION DU CENTRE DE GRAVITE DE MAI2
  130. sx=0.
  131. sy=0.
  132. sz=0.
  133. isd=0
  134. Do 2875 ia=1,nbpts
  135. If (icpr2(ia).ne.0) then
  136. sx=sx+xcoor((ia-1)*(idim+1)+1)
  137. sy=sy+xcoor((ia-1)*(idim+1)+2)
  138. sz=sz+xcoor((ia-1)*(idim+1)+3)
  139. isd = isd+1
  140. Endif
  141. 2875 continue
  142.  
  143. C CREATION DE LA MATRICE DE CHANGEMENT DE BASE
  144. v(1)=(sx/isd)-w(1)
  145. v(2)=(sy/isd)-w(2)
  146. v(3)=(sz/isd)-w(3)
  147. If ((abs(v(1))).lt.(1.E-3) .and. (abs(v(2))).lt.(1.E-3) .and.
  148. $ (abs(v(3))).lt.(1.E-3)) then
  149. v(1)=v(1)+0.1
  150. v(2)=v(2)+0.1
  151. v(3)=v(3)+0.1
  152. Endif
  153. segini mcoor3
  154. segini mcoor1
  155. segini mcoor4
  156. segini mcoor2
  157. r1=sqrt(v(1)*v(1)+v(2)*v(2)+v(3)*v(3))
  158. v(1) = V(1) / r1
  159. v(2) = v(2) / r1
  160. v(3) = v(3) / r1
  161. r2=sqrt(v(1)*v(1)+v(2)*v(2))
  162. if ( r2 . gt. 1.e-9) then
  163. tabmat(1,1)=-v(2)/r2
  164. tabmat(1,2)=v(1)/r2
  165. tabmat(1,3)=0.
  166. else
  167. tabmat(1,1)=1.
  168. tabmat(1,2)=0.
  169. tabmat(1,3)=0.
  170. endif
  171. tabmat(2,1)=-v(3)*tabmat(1,2)
  172. tabmat(2,2)=v(3)*tabmat(1,1)
  173. tabmat(2,3)=v(1)*tabmat(1,2)-v(2)*tabmat(1,1)
  174. tabmat(3,1)=v(1)
  175. tabmat(3,2)=v(2)
  176. tabmat(3,3)=v(3)
  177. * pas besoin de diviser par le determinant car egal à 1.
  178. INVMA(1,1)= TABMAT(2,2)*TABMAT(3,3)-TABMAT(3,2)*TABMAT(2,3)
  179. INVMA(1,2)=-TABMAT(1,2)*TABMAT(3,3)+TABMAT(3,2)*TABMAT(1,3)
  180. INVMA(1,3)= TABMAT(1,2)*TABMAT(2,3)-TABMAT(2,2)*TABMAT(1,3)
  181. INVMA(2,1)=-TABMAT(2,1)*TABMAT(3,3)+TABMAT(3,1)*TABMAT(2,3)
  182. INVMA(2,2)= TABMAT(1,1)*TABMAT(3,3)-TABMAT(3,1)*TABMAT(1,3)
  183. INVMA(2,3)=-TABMAT(1,1)*TABMAT(2,3)+TABMAT(2,1)*TABMAT(1,3)
  184. INVMA(3,1)= TABMAT(2,1)*TABMAT(3,2)-TABMAT(2,2)*TABMAT(3,1)
  185. INVMA(3,2)=-TABMAT(1,1)*TABMAT(3,2)+TABMAT(3,1)*TABMAT(1,2)
  186. INVMA(3,3)= TABMAT(1,1)*TABMAT(2,2)-TABMAT(2,1)*TABMAT(1,2)
  187. C CALCUL DES COORDONNEES DANS LE NOUVEAU REPERE POUR MELEME
  188. nbpts=xcoor(/1)/<ópaj style="color: #009900;">(idim+1)
  189. nbelem=num(/2)
  190. nbnn=num(/1)
  191. segini crit
  192. Do 31 i1=1,nbelem
  193. Do 38 k=1,nbnn
  194. ia=num(k,i1)
  195. xia1=xcoor((ia-1)*(idim+1)+1)-w(1)
  196. xia2=xcoor((ia-1)*(idim+1)+2)-w(2)
  197. xia3=xcoor((ia-1)*(idim+1)+3)-w(3)
  198. Do 32 j=1,3
  199. xcoor3((ia-1)*(idim+1)+j)=((xia1*tabmat(j,1))+
  200. $(xia2*tabmat(j,2))+(xia3*tabmat(j,3)))
  201. 32 continue
  202. 38 continue
  203. 31 continue
  204.  
  205. C CALCUL DE LA TAILLE LA PLUS GRANDE DES ELEMENTS SUIVANT OZ
  206. zlong=-xgrand
  207. Do 1534 i1=1,nbelem
  208. zmin=xgrand
  209. zmax=-xgrand
  210. Do 1233 i2=1,nbnn
  211. ib=num(i2,i1)
  212. ia=(ib-1)*4
  213. If (xcoor3(ia+3).lt.zmin) zmin=xcoor3(ia+3)
  214. If (xcoor3(ia+3).gt.zmax) zmax=xcoor3(ia+3)
  215. 1233 continue
  216. zlong=max(zlong,zmax-zmin)
  217. 1534 continue
  218. If (zlong.lt.(5.E-2)) zlong=(5.E-2)
  219.  
  220. C DETERMINATION DES ELEMENTS TOUCHANT LA TRANCHE CRITIQUE
  221. Do 317 i1=1,nbelem
  222. ime=0
  223. Do 387 k=1,nbnn
  224. ia=num(k,i1)
  225. q(1)=(xcoor3((ia-1)*(idim+1)+1))
  226. q(2)=(xcoor3((ia-1)*(idim+1)+2))
  227. q(3)=(xcoor3((ia-1)*(idim+1)+3))
  228. If (abs(xcoor3((ia-1)*(idim+1)+3)).gt.(0.55*zlong)) then
  229. xcoor1((ia-1)*(idim+1)+1)=(zchoix/(abs(q(3))))*q(1)
  230. xcoor1((ia-1)*(idim+1)+2)=(zchoix/(abs(q(3))))*q(2)
  231. xcoor1((ia-1)*(idim+1)+3)=zchoix
  232. else
  233. ime=ime+1
  234. Endif
  235. 387 continue
  236. If (ime.ne.0) then
  237. tabcrit(1)=tabcrit(1)+1
  238. tabcrit(tabcrit(1)+1)=i1
  239. elecri(i1)=1
  240. Endif
  241. 317 continue
  242.  
  243. C CHANGEMENT DE REPERE POUR IPT1
  244. Do 33 ia=1,nbpts
  245. If (icpr2(ia).eq.1) then
  246. xia1=xcoor((ia-1)*(idim+1)+1)-w(1)
  247. xia2=xcoor((ia-1)*(idim+1)+2)-w(2)
  248. xia3=xcoor((ia-1)*(idim+1)+3)-w(3)
  249. Do 34 j=1,3
  250. xcoor4((ia-1)*(idim+1)+j)=((xia1*tabmat(j,1))+
  251. $(xia2*tabmat(j,2))+(xia3*tabmat(j,3)))
  252. 34 continue
  253. If (abs(xcoor4((ia-1)*(idim+1)+3)).lt.(0.55*zlong)) then
  254. icpr2(ia)=3
  255. else
  256. q(1)=(xcoor4((ia-1)*(idim+1)+1))
  257. q(2)=(xcoor4((ia-1)*(idim+1)+2))
  258. q(3)=(xcoor4((ia-1)*(idim+1)+3))
  259. xcoor2((ia-1)*(idim+1)+1)=(zchoix/(abs(q(3))))*q(1)
  260. xcoor2((ia-1)*(idim+1)+2)=(zchoix/(abs(q(3))))*q(2)
  261. xcoor2((ia-1)*(idim+1)+3)=zchoix
  262. Endif
  263. Endif
  264. If (icpr2(ia).eq.2) then
  265. xia1=xcoor((ia-1)*(idim+1)+1)-w(1)
  266. xia2=xcoor((ia-1)*(idim+1)+2)-w(2)
  267. xia3=xcoor((ia-1)*(idim+1)+3)-w(3)
  268. Do 35 j=1,3
  269. xcoor2((ia-1)*(idim+1)+j)=((xia1*tabmat(j,1))+
  270. $(xia2*tabmat(j,2))+(xia3*tabmat(j,3)))
  271. 35 continue
  272. q(1)=(xcoor4((ia-1)*(idim+1)+1))
  273. q(2)=(xcoor4((ia-1)*(idim+1)+2))
  274. q(3)=(xcoor4((ia-1)*(idim+1)+3))
  275. xcoor2((ia-1)*(idim+1)+1)=(zchoix/(abs(q(3))))*q(1)
  276. xcoor2((ia-1)*(idim+1)+2)=(zchoix/(abs(q(3))))*q(2)
  277. xcoor2((ia-1)*(idim+1)+3)=zchoix
  278. Endif
  279. 33 continue
  280.  
  281. C******************** DEBUT DU ZONAGE DE MAI1 **********************
  282.  
  283. C ON CALCULE LA TAILLE MAXI D'UN ELEMENT DANS TOUTES LES DIRECTIONS
  284. C AFIN DE CREER UN ZONAGE DE L'ESPACE. EN MEME TEMPS ON CALCULE
  285. C LA DIMENSION HORS TOUT DU MAILLAGE
  286. C
  287. IDIM1=4
  288. NBEL = NUM(/2)
  289. NBNN=NUM(/1)
  290. SEGINI ISEG1
  291. ILOC=0
  292. XZO=0.D0
  293. YZO=0.D0
  294. ZZO=0.D0
  295. XZA=XGRAND
  296. YZA=XGRAND
  297. ZZA=XGRAND
  298. XTOMI=XGRAND
  299. XTOMA=-XGRAND
  300. YTOMI=XGRAND
  301. YTOMA=-XGRAND
  302. ZTOMI=XGRAND
  303. ZTOMA=-XGRAND
  304. DO 1 I1=1,NBEL
  305. If (elecri(i1).eq.0) then
  306. XMI=XGRAND
  307. YMI=XGRAND
  308. ZMI=XGRAND
  309. YMA=-XGRAND
  310. XMA=-XGRAND
  311. ZMA=-XGRAND
  312. DO 2 I2 = 1,NBNN
  313. IB=NUM(I2,I1)
  314. IA=(IB-1)*IDIM1
  315. IF(XCOOR1(IA+1).LT.XMI) XMI=XCOOR1(IA+1)
  316. IF(XCOOR1(IA+1).GT.XMA) XMA=XCOOR1(IA+1)
  317. IF(XCOOR1(IA+2).LT.YMI) YMI=XCOOR1(IA+2)
  318. IF(XCOOR1(IA+2).GT.YMA) YMA=XCOOR1(IA+2)
  319. 2 CONTINUE
  320. XLIM(1,I1)=XMI
  321. XLIM(2,I1)=XMA
  322. YLIM(1,I1)=YMI
  323. YLIM(2,I1)=YMA
  324. XZO=MAX (XZO,XMA-XMI)
  325. YZO=MAX (YZO,YMA-YMI)
  326. XZA=MIN(XZA,XMA-XMI)
  327. YZA=MIN(YZA,YMA-YMI)
  328. IF(XMI.LT.XTOMI) XTOMI=XMI
  329. IF(XMA.GT.XTOMA) XTOMA=XMA
  330. IF(YMI.LT.YTOMI) YTOMI=YMI
  331. IF(YMA.GT.YTOMA) YTOMA=YMA
  332. Endif
  333. 1 CONTINUE
  334. XPR=MIN(XZO*1.D-2,(XTOMA-XTOMI)/2.D+4)
  335. YPR=MIN(YZO*1.D-2,(YTOMA-YTOMI)/2.D+4)
  336. XZA=XZA*0.97
  337. YZA=YZA*0.97
  338. XTOMI= XTOMI - (XTOMA-XTOMI)/1.D+4
  339. XTOMA= XTOMA + (XTOMA-XTOMI)/1.D+4
  340. YTOMI= YTOMI - (YTOMA-YTOMI)/1.D+4
  341. YTOMA= YTOMA + (YTOMA-YTOMI)/1.D+4
  342. XZA=MIN ( XZA, XTOMA-XTOMI)
  343. YZA=MIN ( YZA, YTOMA-YTOMI)
  344. XZA = max ( xza , (XTOMA-XTOMI) / 50)
  345. YZA = max ( yza , (YTOMA-YTOMI) / 50)
  346. NXZO=(XTOMA-XTOMI)/XZA + 1
  347. NYZO=(YTOMA-YTOMI)/YZA + 1
  348. XZO=XZA
  349. YZO=YZA
  350. NZZO=1
  351. NXDEP=MIN(NXZO,10)
  352. NYDEP=MIN(NYZO,10)
  353. IF(FLOAT(NXZO)*FLOAT(NYZO).GT.10000.) THEN
  354. XY=SQRT(FLOAT(NXZO)*FLOAT(NYZO))/90
  355. NXZO=MAX(INT(NXZO/XY),NXDEP)
  356. NYZO=MAX(INT(NYZO/XY),NYDEP)
  357. IF(FLOAT(NXZO)*FLOAT(NYZO).GT.10000.) THEN
  358. XY=SQRT(FLOAT(NXZO)*FLOAT(NYZO))/60
  359. NXZO=MAX(INT(NXZO/XY),NXDEP)
  360. NYZO=MAX(INT(NYZO/XY),NYDEP)
  361. ENDIF
  362. XZO=(XTOMA-XTOMI)/NXZO
  363. YZO=(YTOMA-YTOMI)/NYZO
  364. NXZO=(XTOMA-XTOMI)/XZO +1
  365. NYZO=(YTOMA-YTOMI)/YZO +1
  366. ENDIF
  367.  
  368. C
  369. C ON VEUT CONSTRUIRE LA LISTE DES ELEMENTS TOUCHANT UNE ZONE
  370. C POUR CELA ON COMMENCE PAR COMPTER COMBIEN D'ELEMENT TOUCHENT
  371. C CHAQUE ZONE ET EN MEME TEMPS ON STOCKE LES ZONES TOUCHEES
  372. C PAR CHAQUE ELEMENT ET LEUR NOMBRE
  373. C
  374.  
  375. NZO=NXZO*NYZO
  376. IF(IIMPI.NE.0)WRITE(IOIMP,FMT='('' NZO NXZO NYZO NZZO ''
  377. $,4I7) ') NZO,NXZO,NYZO,NZZO
  378. NXYZO=NXZO*NYZO
  379. SEGINI ISEG3
  380. SEGINI ISEG4
  381. DO 3 I1=1,NBEL
  382. If (elecri(i1).eq.0) then
  383. NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) +1
  384. NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) +1
  385. NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) +1
  386. NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) +1
  387. DO 201 L1=NIZ1Y,NIZ2Y
  388. DO 201 L2=NIZ1X,NIZ2X
  389. NIZA = L2 + ( L1-1) * NXZO
  390. NUMZO(NIZA) = NUMZO(NIZA) +1
  391. 201 CONTINUE
  392. Endif
  393. 3 CONTINUE
  394. C
  395. C CONSTRUCTION DU TABLEAU D'ADRESSAGE DU TABLEAU DONNANT LES
  396. C ELEMENTS CONCERNEES PAR UNE ZONE
  397. C
  398. ILON=0
  399. NIZO(1)=1
  400. DO 202 L1=1,NZO
  401. NIZO(L1+1)=NIZO(L1)+NUMZO(L1)
  402. ILON=ILON+ NUMZO(L1)
  403. 202 CONTINUE
  404. 110 FORMAT(16I5)
  405. SEGINI ISEG5
  406. DO 5 I1=1,NBEL
  407. If (elecri(i1).eq.0) then
  408. NIZ1X=INT((XLIM(1,I1)-XTOMI-XPR)/XZO) +1
  409. NIZ1Y=INT((YLIM(1,I1)-YTOMI-YPR)/YZO) +1
  410. NIZ2X=INT((XLIM(2,I1)-XTOMI+XPR)/XZO) +1
  411. NIZ2Y=INT((YLIM(2,I1)-YTOMI+YPR)/YZO) +1
  412. DO 203 L1=NIZ1Y,NIZ2Y
  413. DO 203 L2=NIZ1X,NIZ2X
  414. NIZA = L2 + ( L1-1) * NXZO
  415. IAD=NIZO(NIZA)+IDEJ(NIZA)
  416. NNMEL(IAD)=I1
  417. IDEJ(NIZA)=IDEJ(NIZA)+1
  418. 203 CONTINUE
  419. Endif
  420. 5 CONTINUE
  421. C *********************** FIN DU ZONAGE DE MAI1 **************************
  422. prec = 1.E-5
  423. prec2 = -prec
  424. prec3 = -5.E-3
  425. C ******RECHERCHE DU NOMBRE D'ELEMENT MAX CONTENU DANS UNE ZONE***********
  426. la=0
  427. ls=0
  428. Do 420 i=1,nzo
  429. ls=ls+numzo(i)
  430. If (numzo(i).gt.la) then
  431. la=numzo(i)
  432. Endif
  433. 420 continue
  434. ionmax=la
  435.  
  436. C ****************** TRAITEMENT SUR LES POINTS DE MAI2 *******************
  437. nbelem=num(/2)
  438. nbnn=num(/1)
  439. nbpts=xcoor(/1)/(idim+1)
  440. nbpts2=nbpts
  441. segini mcor,mcorre,iseg6
  442. Do 405 ia=1,nbpts2
  443. If (icpr2(ia).eq.2) then
  444. icorre(ia)=ia
  445. Endif
  446. If (icpr2(ia).eq.1) then
  447. xpu(1)=xcoor2((ia-1)*(idim+1)+1)
  448. xpu(2)=xcoor2((ia-1)*(idim+1)+2)
  449. If (tabcrit(1).eq.0) then
  450. if (((xpu(1)-xtomi).lt.prec.) .or.
  451. $ ((xpu(2)-ytomi).lt.prec.) .or.
  452. $ ((xpu(1)-xtoma).gt.prec2.) .or.
  453. $ ((xpu(2)-ytoma).gt.prec2.)) then
  454. call erreur(21)
  455. return
  456. endif
  457. Endif
  458. C RECHERCHE DU NUMERO DE LA ZONE CORRESPONDANTE
  459. k2=int(((xcoor2((ia-1)*(idim+1)+1))-xtomi-xpr)/xzo)+1
  460. k1=int(((xcoor2((ia-1)*(idim+1)+2))-ytomi-ypr)/yzo)+1
  461. If ((k2.gt.0.) .and. (k2.lt.(nxzo+1)) .and.
  462. $ (k1.gt.0.) .and. (k1.lt.(nyzo+1))) then
  463. niza=k2+(k1-1)*nxzo
  464. C INVENTAIRE DES ELEMENTS CONCERNES PAR LE POINT IA
  465. iad=nizo(niza)
  466. Do 406 i=1,numzo(niza)
  467. i1=nnmel(iad)
  468. C CALCUL DES COORDONNEES BARYCENTRIQUES DE IA DANS L'ELEMENT I1
  469. j1=num(1,i1)
  470. j2=num(2,i1)
  471. j3=num(3,i1)
  472. j1idim=(j1-1)*(idim+1)
  473. j2idim=(j2-1)*(idim+1)
  474. j3idim=(j3-1)*(idim+1)
  475. Do 408 l=1,2
  476. p(l+1)=xpu(l)
  477. am(l+1,1)=xcoor1(j1idim+l)
  478. am(l+1,2)=xcoor1(j2idim+l)
  479. am(l+1,3)=xcoor1(j3idim+l)
  480. 408 continue
  481. x1=am(2,1)
  482. x2=am(2,2)
  483. x3=am(2,3)
  484. y1=am(3,1)
  485. y2=am(3,2)
  486. y3=am(3,3)
  487. x=p(2)
  488. y=p(3)
  489. detam=x1*y2+x2*y3+x3*y1-y1*x2-y2*x3-y3*x1
  490. a(1)=x2*y3-y2*x3
  491. a(2)=x3*y1-x1*y3
  492. a(3)=x1*y2-x2*y1
  493. b(1)=y2-y3
  494. b(2)=y3-y1
  495. b(3)=y1-y2
  496. c(1)=x3-x2
  497. c(2)=x1-x3
  498. c(3)=x2-x1
  499. Do 409 k=1,nbnn
  500. al(k)=(a(k)+b(k)*x+c(k)*y)/detam
  501. 409 continue
  502. C TEST DE POSITION SUR LES COORDONNEES BARYCENTRIQUES
  503. If (al(1).gt.prec3 . and. al(2).gt.prec3 . and.
  504. $ al(3).gt.prec3) then
  505. C LE POINT EST A L'INTERIEUR DE L'ELEMENT I1
  506. icor(1)=icor(1)+1
  507. icor(icor(1)+1)=i1
  508. Endif
  509. iad=iad+1
  510.  
  511. 406 continue
  512. Endif
  513. C PASSAGE TRANCHE CRITIQUE POUR UN POINT DU ZONAGE
  514. If (k2.lt.2 .or. k2.gt.(nxzo-1)) then
  515. Do 1462 iz=1,tabcrit(1)
  516. i1=tabcrit(iz+1)
  517. m1=num(1,i1)
  518. m2=num(2,i1)
  519. m3=num(3,i1)
  520. C CALCUL DES COORDONNEES DU PROJETE DE IA SUR L'ELEMENT ICOR(2)
  521. xm1=xcoor3((m1-1)*(idim+1)+1)
  522. ym1=xcoor3((m1-1)*(idim+1)+2)
  523. zm1=xcoor3((m1-1)*(idim+1)+3)
  524. xm2=xcoor3((m2-1)*(idim+1)+1)
  525. ym2=xcoor3((m2-1)*(idim+1)+2)
  526. zm2=xcoor3((m2-1)*(idim+1)+3)
  527. xm3=xcoor3((m3-1)*(idim+1)+1)
  528. ym3=xcoor3((m3-1)*(idim+1)+2)
  529. zm3=xcoor3((m3-1)*(idim+1)+3)
  530. q(1)=(xcoor4((ia-1)*(idim+1)+1))
  531. q(2)=(xcoor4((ia-1)*(idim+1)+2))
  532. q(3)=(xcoor4((ia-1)*(idim+1)+3))
  533. C Calcul du zib
  534. o1=(ym1-ym2)*(zm1-zm3)-(ym1-ym3)*(zm1-zm2)
  535. o2=(zm1-zm2)*(xm1-xm3)-(zm1-zm3)*(xm1-xm2)
  536. o3=(xm1-xm2)*(ym1-ym3)-(xm1-xm3)*(ym1-ym2)
  537. o4=-(o1*xm1+o2*ym1+o3*zm1)
  538. zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3)
  539. xnew=q(1)*zk
  540. ynew=q(2)*zk
  541. znew=q(3)*zk
  542. C PASSAGE DANS UN REPERE DONT L'AXE Z EST ORTHO AU PLAN DE I1
  543. r1=sqrt(o1*o1+o2*o2+o3*o3)
  544. if( r1 . EQ. 0. ) then
  545. call erreur ( 21 )
  546. endif
  547. o1 = o1/r1
  548. o2 = o2/r1
  549. o3 = o3/r1
  550. r2=sqrt(o1*o1+o2*o2)
  551. if ( r2 . gt. 1.e-5) then
  552. tabmat(1,1)=-o2/r2
  553. tabmat(1,2)=o1/r2
  554. tabmat(1,3)=0.
  555. else
  556. tabmat(1,1)=1.
  557. tabmat(1,2)=0.
  558. tabmat(1,3)=0.
  559. endif
  560. tabmat(2,1)=-o3*tabmat(1,2)
  561. tabmat(2,2)=o3*tabmat(1,1)
  562. tabmat(2,3)=o1*tabmat(1,2)-o2*tabmat(1,1)
  563. tabmat(3,1)=o1
  564. tabmat(3,2)=o2
  565. tabmat(3,3)=o3
  566. x=((xnew*tabmat(1,1))+(ynew*tabmat(1,2))+(znew*tabmat(1,3)))
  567. y=((xnew*tabmat(2,1))+(ynew*tabmat(2,2))+(znew*tabmat(3,3)))
  568. x1=((xm1*tabmat(1,1))+(ym1*tabmat(1,2))+(zm1*tabmat(1,3)))
  569. y1=((xm1*tabmat(2,1))+(ym1*tabmat(2,2))+(zm1*tabmat(3,3)))
  570. x2=((xm2*tabmat(1,1))+(ym2*tabmat(1,2))+(zm2*tabmat(1,3)))
  571. y2=((xm2*tabmat(2,1))+(ym2*tabmat(2,2))+(zm2*tabmat(3,3)))
  572. x3=((xm3*tabmat(1,1))+(ym3*tabmat(1,2))+(zm3*tabmat(1,3)))
  573. y3=((xm3*tabmat(2,1))+(ym3*tabmat(2,2))+(zm3*tabmat(3,3)))
  574. C CALCUL DES COORDONNEES BARYCENTRIQUES
  575. detam=x1*y2+x2*y3+x3*y1-y1*x2-y2*x3-y3*x1
  576. a(1)=x2*y3-y2*x3
  577. a(2)=x3*y1-x1*y3
  578. a(3)=x1*y2-x2*y1
  579. b(1)=y2-y3
  580. b(2)=y3-y1
  581. b(3)=y1-y2
  582. c(1)=x3-x2
  583. c(2)=x1-x3
  584. c(3)=x2-x1
  585. Do 798 k=1,nbnn
  586. al(k)=(a(k)+b(k)*x+c(k)*y)/detam
  587. 798 continue
  588. C TEST DE POSITION SUR LES COORDONNEES BARYCENTRIQUES
  589. If (al(1).gt.prec3 . and. al(2).gt.prec3 . and.
  590. $ al(3).gt.prec3) then
  591. C LE POINT EST A L'INTERIEUR DE L'ELEMENT I1
  592. icor(1)=icor(1)+1
  593. icor(icor(1)+1)=i1
  594. Endif
  595. 1462 continue
  596. Endif
  597.  
  598. If (icor(1).eq.0) then
  599. interr(1) = ia
  600. call erreur(782)
  601. return
  602. Endif
  603.  
  604. Endif
  605. C FIN DU IF ICPR2(IA).EQ.1
  606.  
  607. If (icpr2(ia).eq.3) then
  608.  
  609. C PASSAGE SUR LES ELEMENTS DE LA TRANCHE CRITIQUE
  610. Do 1463 iz=1,tabcrit(1)
  611. i1=tabcrit(iz+1)
  612. m1=num(1,i1)
  613. m2=num(2,i1)
  614. m3=num(3,i1)
  615. C CALCUL DES COORDONNEES DU PROJETE DE IA SUR L'ELEMENT ICOR(2)
  616. xm1=xcoor3((m1-1)*(idim+1)+1)
  617. ym1=xcoor3((m1-1)*(idim+1)+2)
  618. zm1=xcoor3((m1-1)*(idim+1)+3)
  619. xm2=xcoor3((m2-1)*(idim+1)+1)
  620. ym2=xcoor3((m2-1)*(idim+1)+2)
  621. zm2=xcoor3((m2-1)*(idim+1)+3)
  622. xm3=xcoor3((m3-1)*(idim+1)+1)
  623. ym3=xcoor3((m3-1)*(idim+1)+2)
  624. zm3=xcoor3((m3-1)*(idim+1)+3)
  625. q(1)=(xcoor4((ia-1)*(idim+1)+1))
  626. q(2)=(xcoor4((ia-1)*(idim+1)+2))
  627. q(3)=(xcoor4((ia-1)*(idim+1)+3))
  628. C Calcul du zib
  629. o1=(ym1-ym2)*(zm1-zm3)-(ym1-ym3)*(zm1-zm2)
  630. o2=(zm1-zm2)*(xm1-xm3)-(zm1-zm3)*(xm1-xm2)
  631. o3=(xm1-xm2)*(ym1-ym3)-(xm1-xm3)*(ym1-ym2)
  632. o4=-(o1*xm1+o2*ym1+o3*zm1)
  633. zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3)
  634. xnew=q(1)*zk
  635. ynew=q(2)*zk
  636. znew=q(3)*zk
  637. C PASSAGE DANS UN REPERE DONT L'AXE Z EST ORTHO AU PLAN DE I1
  638. r1=sqrt(o1*o1+o2*o2+o3*o3)
  639. if( r1 . EQ. 0. ) then
  640. call erreur ( 21 )
  641. endif
  642. o1 = o1/r1
  643. o2 = o2/r1
  644. o3 = o3/r1
  645. r2=sqrt(o1*o1+o2*o2)
  646. if ( r2 . gt. 1.e-5) then
  647. tabmat(1,1)=-o2/r2
  648. tabmat(1,2)=o1/r2
  649. tabmat(1,3)=0.
  650. else
  651. tabmat(1,1)=1.
  652. tabmat(1,2)=0.
  653. tabmat(1,3)=0.
  654. endif
  655. tabmat(2,1)=-o3*tabmat(1,2)
  656. tabmat(2,2)=o3*tabmat(1,1)
  657. tabmat(2,3)=o1*tabmat(1,2)-o2*tabmat(1,1)
  658. tabmat(3,1)=o1
  659. tabmat(3,2)=o2
  660. tabmat(3,3)=o3
  661. x=((xnew*tabmat(1,1))+(ynew*tabmat(1,2))+(znew*tabmat(1,3)))
  662. y=((xnew*tabmat(2,1))+(ynew*tabmat(2,2))+(znew*tabmat(3,3)))
  663. x1=((xm1*tabmat(1,1))+(ym1*tabmat(1,2))+(zm1*tabmat(1,3)))
  664. y1=((xm1*tabmat(2,1))+(ym1*tabmat(2,2))+(zm1*tabmat(3,3)))
  665. x2=((xm2*tabmat(1,1))+(ym2*tabmat(1,2))+(zm2*tabmat(1,3)))
  666. y2=((xm2*tabmat(2,1))+(ym2*tabmat(2,2))+(zm2*tabmat(3,3)))
  667. x3=((xm3*tabmat(1,1))+(ym3*tabmat(1,2))+(zm3*tabmat(1,3)))
  668. y3=((xm3*tabmat(2,1))+(ym3*tabmat(2,2))+(zm3*tabmat(3,3)))
  669. C CALCUL DES COORDONNEES BARYCENTRIQUES
  670. detam=x1*y2+x2*y3+x3*y1-y1*x2-y2*x3-y3*x1
  671. a(1)=x2*y3-y2*x3
  672. a(2)=x3*y1-x1*y3
  673. a(3)=x1*y2-x2*y1
  674. b(1)=y2-y3
  675. b(2)=y3-y1
  676. b(3)=y1-y2
  677. c(1)=x3-x2
  678. c(2)=x1-x3
  679. c(3)=x2-x1
  680. Do 799 k=1,nbnn
  681. al(k)=(a(k)+b(k)*x+c(k)*y)/detam
  682. 799 continue
  683. C TEST DE POSITION SUR LES COORDONNEES BARYCENTRIQUES
  684. If (al(1).gt.prec3 . and. al(2).gt.prec3 . and.
  685. $al(3).gt.prec3) then
  686. C LE POINT EST A L'INTERIEUR DE L'ELEMENT I1
  687. icor(1)=icor(1)+1
  688. icor(icor(1)+1)=i1
  689. Endif
  690. 1463 continue
  691. If (icor(1).eq.0) then
  692. interr(1)=ia
  693. call erreur(782)
  694. return
  695. Endif
  696. Endif
  697. C FIN DE L'INVENTAIRE
  698.  
  699. C CALCUL DES PROJETES DE MAI2 SUR MAI1
  700. If ((icpr2(ia).eq.1) .or. (icpr2(ia).eq.3)) then
  701.  
  702. dmin=1.e30
  703. nbpts=nbpts+1
  704. segadj mcoord
  705. n=0
  706. Do 500 m=1,icor(1)
  707. m1=num(1,icor(m+1))
  708. m2=num(2,icor(m+1))
  709. m3=num(3,icor(m+1))
  710. C Calcul des coordonnees du projete
  711. xm1=xcoor3((m1-1)*(idim+1)+1)
  712. ym1=xcoor3((m1-1)*(idim+1)+2)
  713. zm1=xcoor3((m1-1)*(idim+1)+3)
  714. xm2=xcoor3((m2-1)*(idim+1)+1)
  715. ym2=xcoor3((m2-1)*(idim+1)+2)
  716. zm2=xcoor3((m2-1)*(idim+1)+3)
  717. xm3=xcoor3((m3-1)*(idim+1)+1)
  718. ym3=xcoor3((m3-1)*(idim+1)+2)
  719. zm3=xcoor3((m3-1)*(idim+1)+3)
  720. q(1)=(xcoor4((ia-1)*(idim+1)+1))
  721. q(2)=(xcoor4((ia-1)*(idim+1)+2))
  722. q(3)=(xcoor4((ia-1)*(idim+1)+3))
  723. C Calcul du zib
  724. o1=(ym1-ym2)*(zm1-zm3)-(ym1-ym3)*(zm1-zm2)
  725. o2=(zm1-zm2)*(xm1-xm3)-(zm1-zm3)*(xm1-xm2)
  726. o3=(xm1-xm2)*(ym1-ym3)-(xm1-xm3)*(ym1-ym2)
  727. o4=-(o1*xm1+o2*ym1+o3*zm1)
  728. zk=-o4/(q(1)*o1+q(2)*o2+q(3)*o3)
  729. xnew=q(1)*zk
  730. ynew=q(2)*zk
  731. znew=q(3)*zk
  732. delta1=sqrt(xnew*xnew+ynew*ynew+znew*znew)
  733. delta2=sqrt(q(1)*q(1)+q(2)* q(2)+q(3)*q(3))
  734. pds=(xnew*q(1))+(ynew*q(2))+(znew*q(3))
  735. If (pds.gt.0.) then
  736. If ((delta1-delta2).gt.0.) then
  737. If ((delta1-delta2).lt.dmin) then
  738. icorre(ia) = nbpts
  739. dmin=(delta1-delta2)
  740. xcoor((nbpts-1)*(idim+1)+1)=xnew
  741. xcoor((nbpts-1)*(idim+1)+2)=ynew
  742. xcoor((nbpts-1)*(idim+1)+3)=znew
  743. Endif
  744. n=n+1
  745. Endif
  746. Endif
  747. 500 continue
  748. If (n.eq.0) then
  749. call erreur(21)
  750. return
  751. Endif
  752.  
  753. C *********************RETOUR A L'ANCIENNE BASE*********************
  754. ib=nbpts
  755. xa1=xcoor((ib-1)*(idim+1)+1)
  756. xa2=xcoor((ib-1)*(idim+1)+2)
  757. xa3=xcoor((ib-1)*(idim+1)+3)
  758. do 3200 j=1,3
  759. xcoor((ib-1)*(idim+1)+j)=(xa1*invma(j,1))+
  760. $(xa2*invma(j,2))+(xa3*invma(j,3))+w(J)
  761. 3200 continue
  762. icor(1)=0
  763. n=0
  764. Endif
  765. 405 continue
  766. segsup iseg1,iseg3,iseg4,iseg5,iseg6,mcoor1,mcoor2
  767. segsup icp1,icp2,mcoor3,mcoor4
  768. C CREATION DU MAILLAGE PROJETE MAI3
  769. segini,ipt3=ipt1
  770. ipt4 = IPT3
  771. do 2456 IU = 1, max( 1,IPT3.LISOUS(/1))
  772. If( IPT3.LISOUS(/1).NE.0) THEN
  773. IPT5=IPT3.LISOUS(IU)
  774. SEGINI,IPT4=IPT5
  775. NBREF = 0
  776. NBSOUS=0
  777. SEGADJ IPT4
  778. ENDIF
  779. nbel1=ipt4.num(/2)
  780. nbnn1=ipt4.num(/1)
  781. Do 600 i=1,nbel1
  782. Do 601 j=1,nbnn1
  783. ia=ipt4.num(j,i)
  784. ib=icorre(ia)
  785. ipt4.num(j,i)=ib
  786. 601 continue
  787. 600 continue
  788. 2456 continue
  789. do 2457 IU = 1, max( 0,IPT3.LISREF(/1))
  790. If( IPT3.LISREF(/1).NE.0) THEN
  791. IPT5 = IPT3.LISREF(IU)
  792. SEGINI,IPT4=IPT5
  793. ENDIF
  794. nbel1=ipt4.num(/2)
  795. nbnn1=ipt4.num(/1)
  796. Do 602 i=1,nbel1
  797. Do 603 j=1,nbnn1
  798. ia=ipt4.num(j,i)
  799. ib=icorre(ia)
  800. ipt4.num(j,i)=ib
  801. 603 continue
  802. 602 continue
  803. SEGDES IPT4
  804. 2457 continue
  805. call ecrobj('MAILLAGE',ipt3)
  806. segdes ipt3
  807. segsup mcor,mcorre,matric
  808. Return
  809. End
  810.  
  811.  
  812.  
  813.  
  814.  
  815.  
  816.  
  817.  
  818.  
  819.  

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