Télécharger facaxi.eso

Retour à la liste

Numérotation des lignes :

  1. C FACAXI SOURCE CB215821 16/04/21 21:16:44 8920
  2. SUBROUTINE FACAXI (MYMOD,INFOEL,NPAX,NGAX,KACHE,INOR,ICHFAC
  3. & ,EXTINC)
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8 (A-H,O-Z)
  6. C----------------------------------------------------------------------
  7. C SP appele par FFORME
  8. C
  9. C Calcul des facteurs de forme en axisymetrique
  10. C entée:
  11. C MYMOD : pointeur de l'objet modèle
  12. C INFOEL : utile pour les coques ou quadratiques
  13. C NPAX : nombre de points d integration (disc.reguliere)
  14. C NGAX : nombre de points de Gauss
  15. C KACHE : 0 si option convexe, sinon option cache
  16. C INOR : 1 si normalisation, 0 sinon
  17. C EXTINC : coefficient d'extinction (si cavite absorbante)
  18. C sortie:
  19. C ICHFAC : pointeur sur l'objet MCHAML resultat
  20. C
  21. C----------------------------------------------------------------------
  22. C traitement des coques par dedoublement des elements a partir
  23. C de la normale
  24. C ->
  25. C A (inverse) = A - e*n (e=1e-3)
  26. C cas des boucles 1 sur k1 et 2 sur k2
  27. C mais pas de la boucle 3 obstructeurs
  28. C
  29. C bcl face k1
  30. C ** face k1 **
  31. C bcl face k2
  32. C .. face k2 ..
  33. C bcl 3 obstructeurs
  34. C .. si coq: inverse face k2 ..
  35. C les obstructeurs sont les memes que pour k2
  36. C fin bcl face k2
  37. C
  38. C ** si coq : inverse face k1 **
  39. C bcl face k2
  40. C .. face k2 ..
  41. C bcl 3 obstructeurs
  42. C .. si coq: inverse face k2 ..
  43. C les obstructeurs sont les memes que pour k2
  44. C fin bcl face k2
  45. C fin bcl face k1
  46. C----------------------------------------------------------------------
  47. C
  48. LOGICAL ICOQ
  49. C
  50. -INC CCREEL
  51. -INC SMELEME
  52. -INC SMMODEL
  53. -INC CCOPTIO
  54. -INC SMCOORD
  55. POINTEUR MYMOD.MMODEL
  56.  
  57. C----------------------------------------------------------------------
  58. SEGMENT , INFOEL
  59. LOGICAL KCOQ(N1),KQUAD(N1)
  60. ENDSEGMENT
  61. C----------------------------------------------------------------------
  62. C FACTEURS DE FORME
  63. C NNBEL1 = NOMBRE DE LIGNES + 1
  64. C NBEL2 = NOMBRE DE COLONNES
  65. C LFACT(NNBEL1) POINTE SUR LE TABLEAU DES SURFACES
  66. C
  67. SEGMENT IFACFO
  68. INTEGER LFACT(NNBEL1)
  69. ENDSEGMENT
  70. SEGMENT LFAC
  71. REAL*8 FACT(NBEL2)
  72. ENDSEGMENT
  73. POINTEUR PSUR.LFAC, PLIG.LFAC
  74. POINTEUR PCOL.LFAC
  75. C----------------------------------------------------------------------
  76. C coordonnees des obstructeurs
  77. SEGMENT SFOBS
  78. REAL*8 OBS(2,NTOBS)
  79. ENDSEGMENT
  80. C----------------------------------------------------------------------
  81. SEGMENT STRAV
  82. REAL*8 A1(NA,NA),A2(NA,NA),A3(NA,NA),AA2(NA,NA)
  83. REAL*8 U1(NA1),U2(NA1),UU2(NA1)
  84. ENDSEGMENT
  85. C----------------------------------------------------------------------
  86. C
  87. EPS = 1D-5
  88.  
  89. KIMP = IIMPI
  90. NES = IDIM
  91.  
  92. C... critere de dedoublement des coques
  93. ECOQ=1.D-3
  94. IF (INFOEL.EQ.0) THEN
  95. ICOQ = .FALSE.
  96. ELSE
  97. ICOQ = .TRUE.
  98. SEGACT INFOEL
  99. ENDIF
  100. C... quadratique
  101. NSPA1=1
  102. NSPA2=1
  103. NSPA3=1
  104. NS=2
  105.  
  106. RAD = 0
  107. C
  108. C Calcul du nombre de faces NFACE et du rayon RAD
  109. SEGACT MYMOD
  110. NTYP = MYMOD.KMODEL(/1)
  111. NFACE = 0
  112. DO 10 ITYP=1,NTYP
  113. IMODEL = MYMOD.KMODEL(ITYP)
  114. SEGACT IMODEL
  115. IPT1 = IMAMOD
  116. SEGDES IMODEL
  117. SEGACT IPT1
  118. NEL = IPT1.NUM(/2)
  119. NSGEO = IPT1.NUM(/1)
  120. C Recherche du max sur Ox
  121. DO 5 IEL = 1,NEL
  122. DO 6 IPT = 1,NSGEO
  123. IREF = (IDIM+1)*(IPT1.NUM(IPT,IEL)-1)
  124. VALX = XCOOR(IREF+1)
  125. IF (VALX.GT.RAD) RAD = VALX
  126. 6 CONTINUE
  127. 5 CONTINUE
  128. IF (ICOQ.AND.KCOQ(ITYP)) THEN
  129. NFACE = NFACE + 2 * NEL
  130. ELSE
  131. NFACE = NFACE + NEL
  132. ENDIF
  133.  
  134. 10 CONTINUE
  135.  
  136. C
  137.  
  138. IF (KIMP.GE.3) THEN
  139. WRITE( 6,*) ' DIMENSIONNEMENT : RAD =',RAD
  140. WRITE( 6,*) ' NOMBRE TOTAL DE FACES ',NFACE
  141. ENDIF
  142. C
  143. C>>> INITIALISATION SFOBS et STRAV
  144. C -----------------------------
  145. IF(KACHE.NE.0) THEN
  146. NTOBS = 2*NFACE
  147. SEGINI SFOBS
  148. ENDIF
  149. NA=2
  150. NA1=3
  151. SEGINI STRAV
  152.  
  153. C
  154. C>>> INITIALISATION OBJET FACFOR
  155. C ---------------------------
  156. C
  157. NNBEL1=NFACE+1
  158. NBEL2=NFACE
  159. SEGINI IFACFO
  160. DO 900 LS=1,NNBEL1
  161. SEGINI PLIG
  162. LFACT(LS)=PLIG
  163. SEGACT PLIG*MOD
  164. 900 CONTINUE
  165. PSUR = LFACT(NNBEL1)
  166. SEGACT PSUR*MOD
  167.  
  168. C -------------------------------------------------------------
  169. C>>> CALCUL
  170.  
  171. C
  172. C>> BOUCLE FACE 1
  173. C
  174.  
  175. NELT1= 0
  176. DO 100 ITYP = 1,NTYP
  177. IMODEL = MYMOD.KMODEL(ITYP)
  178. SEGACT IMODEL
  179. IPT1 = IMAMOD
  180. SEGDES IMODEL
  181. IF(KIMP.GE.4) THEN
  182. WRITE (6,*) 'Maillage :',IPT1
  183. ENDIF
  184.  
  185. NSGEO1 = IPT1.NUM(/1)
  186. NSPA1=1
  187. IF(ICOQ.AND.KQUAD(ITYP)) NSPA1=2
  188. NEL1 = IPT1.NUM(/2)
  189. DO 110 I1 = 1,NEL1
  190.  
  191. IF (ICOQ.AND.KCOQ(ITYP)) THEN
  192. K1 = (2*I1-1) + NELT1
  193. ELSE
  194. K1 = I1+ NELT1
  195. ENDIF
  196.  
  197. C*** traitement de la face K1 ***************************************
  198.  
  199. PLIG=LFACT(K1)
  200. C
  201. DO 111 IS = 1,NSGEO1,NSPA1
  202. LS=IS
  203. IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  204. IREF = (IDIM+1)*(IPT1.NUM(IS,I1)-1)
  205. IF(KIMP.GE.4) THEN
  206. WRITE (6,*) 'Noeud numéro',IPT1.NUM(IS,I1),'ref :'
  207. $ ,IREF
  208. ENDIF
  209. DO 112 K = 1,NES
  210. C** A1(K,IS) = XCOOR(IREF+K)/RAD
  211. A1(K,LS) = XCOOR(IREF+K)/RAD
  212. 112 CONTINUE
  213. 111 CONTINUE
  214. CALL KNORM2(NES,A1,NS,U1,S1)
  215. S1=XPI*S1*ABS(A1(1,1)+A1(1,2))
  216. PSUR.FACT(K1)=S1
  217. ZMIN1=DMIN1(A1(2,1),A1(2,2))
  218. ZMAX1=DMAX1(A1(2,1),A1(2,2))
  219. RMAX1=DMAX1(A1(1,1),A1(1,2))
  220.  
  221.  
  222. C>> BOUCLE FACE 2
  223. C -------------------------------------------------------------
  224.  
  225. NELT2= 0
  226. DO 200 JTYP = 1,NTYP
  227. IMODEL = MYMOD.KMODEL(JTYP)
  228. SEGACT IMODEL
  229. IPT2 = IMAMOD
  230. SEGDES IMODEL
  231.  
  232. NSGEO2 = IPT2.NUM(/1)
  233. NSPA2=1
  234. IF(ICOQ.AND.KQUAD(JTYP)) NSPA2=2
  235. NEL2 = IPT2.NUM(/2)
  236. DO 210 I2 = 1,NEL2
  237.  
  238. IF (ICOQ.AND.KCOQ(JTYP))THEN
  239. K2 = 2*I2-1 + NELT2
  240. ELSE
  241. K2 = I2 + NELT2
  242. ENDIF
  243.  
  244. C... face K2 .....................................................
  245.  
  246. KVU=1
  247. FF=0.D0
  248.  
  249. C.. UTILISATION DE LA RECIPROCITE
  250. C
  251. IF(K1.GT.K2) THEN
  252. S2=PSUR.FACT(K2)
  253. PCOL=LFACT(K2)
  254.  
  255. PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1)
  256.  
  257. ELSE
  258. C.. TEST K1=K2 ET VISIBILITE A PRIORI
  259.  
  260. C------------------------------------------------------------------
  261.  
  262. C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0
  263. C
  264. IF(K1.EQ.K2) THEN
  265. IF(KIMP.GE.4) THEN
  266. WRITE(6,*) ' A1 ',A1(1,1),A1(2,1)
  267. WRITE(6,*) ' A1 ',A1(1,2),A1(2,2)
  268. WRITE(6,*) ' U1 ',U1(1),U1(2)
  269. ENDIF
  270. IF (ABS(U1(1)).LT.EPS) KVU=0
  271. IF (U1(1).GE.EPS) KVU=0
  272.  
  273. IF(KVU.NE.0) THEN
  274. DO 214 K = 1,NES
  275. A2(K,1) = A1(K,1)
  276. A2(K,2) = A1(K,2)
  277. 214 CONTINUE
  278. ENDIF
  279.  
  280. C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE
  281. C UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z
  282. C ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z)
  283.  
  284. ELSE
  285. DO 211 IS = 1,NSGEO2,NSPA2
  286. LS=IS
  287. IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  288. IREF = (IDIM+1)*(IPT2.NUM(IS,I2)-1)
  289. DO 212 K = 1,NES
  290. C** A2(K,IS) = XCOOR(IREF+K)/RAD
  291. A2(K,LS) = XCOOR(IREF+K)/RAD
  292. 212 CONTINUE
  293. 211 CONTINUE
  294. CALL KNORM2(NES,A2,NS,U2,S2)
  295. S2=XPI*S2*ABS(A2(1,1)+A2(1,2))
  296.  
  297. C.. calcul du symetrique
  298. AA2(1,1)= -A2(1,2)
  299. AA2(2,1)= A2(2,2)
  300. AA2(1,2)= -A2(1,1)
  301. AA2(2,2)= A2(2,1)
  302. CALL KNORM2(NES,AA2,NS,UU2,SS2)
  303.  
  304. C.. visibilite a priori
  305. CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU)
  306. CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS)
  307.  
  308. IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS ',K1,K2
  309. $ ,IVU,IVUS
  310. IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0
  311.  
  312. ENDIF
  313. C ---------------------------------------------------------------
  314. C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI
  315. IF(KIMP.GE.4.AND.KVU.NE.0) THEN
  316. WRITE(6,*) ' FACES VISIBLES ',K1,K2
  317. ENDIF
  318.  
  319. C>> CALCUL
  320. C... option convexe
  321. IF(KACHE.EQ.0) THEN
  322. IF(KVU.NE.0) THEN
  323.  
  324. CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP,EXTINC,RAD)
  325. IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE ',K1,K2
  326. $ ,FF/S1
  327.  
  328. ENDIF
  329. C.. option cache
  330. ELSE
  331. IF(KVU.NE.0) THEN
  332.  
  333. C>> RECHERCHE DES OBSTRUCTEURS POTENTIELS.--------------------------
  334.  
  335. NOBS=0
  336.  
  337. C>> BOUCLE FACE 3---------------------------------------------------
  338. NELT3= 0
  339. DO 300 KTYP = 1,NTYP
  340. IMODEL = MYMOD.KMODEL(KTYP)
  341. SEGACT IMODEL
  342. IPT3 = IMAMOD
  343. SEGDES IMODEL
  344.  
  345. NSGEO3 = IPT3.NUM(/1)
  346. NSPA3=1
  347. IF(ICOQ.AND.KQUAD(KTYP)) NSPA3=2
  348. NEL3 = IPT3.NUM(/2)
  349. DO 310 I3 = 1,NEL3
  350. K3 = I3+ NELT3
  351.  
  352. IF(K3.NE.K1.AND.K3.NE.K2) THEN
  353.  
  354. DO 311 IS = 1,NSGEO3,NSPA3
  355. LS=IS
  356. IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1
  357. $ )/2
  358. IREF = (IDIM+1)*(IPT3.NUM(IS,I3)
  359. $ -1)
  360. DO 312 K = 1,NES
  361. C** A3(K,IS) = XCOOR(IREF+K)/RAD
  362. A3(K,LS) = XCOOR(IREF+K)/RAD
  363. 312 CONTINUE
  364. 311 CONTINUE
  365.  
  366. ZMIN2=DMIN1(A2(2,1),A2(2,2))
  367. ZMAX2=DMAX1(A2(2,1),A2(2,2))
  368. RMAX2=DMAX1(A2(1,1),A2(1,2))
  369. ZTMIN=DMIN1(ZMIN1,ZMIN2)
  370. ZTMAX=DMAX1(ZMAX1,ZMAX2)
  371. RMAX =DMAX1(RMAX1,RMAX2)
  372.  
  373. IF(DMAX1(A3(2,1),A3(2,2)).LE.ZTMIN)
  374. $ THEN
  375. KOBS=0
  376. ELSEIF(DMIN1(A3(2,1),A3(2,2)).GE
  377. $ .ZTMAX) THEN
  378. KOBS=0
  379. ELSEIF(DMIN1(A3(1,1),A3(1,2)).GE
  380. $ .RMAX) THEN
  381. KOBS=0
  382. ELSE
  383. KOBS=1
  384. OBS(1,2*NOBS+1)=A3(1,1)
  385. OBS(2,2*NOBS+1)=A3(2,1)
  386. OBS(1,2*NOBS+2)=A3(1,2)
  387. OBS(2,2*NOBS+2)=A3(2,2)
  388. NOBS=NOBS+1
  389.  
  390. ENDIF
  391. ENDIF
  392.  
  393. 310 CONTINUE
  394. C
  395. NELT3=NELT3+NEL3
  396. C
  397. 300 CONTINUE
  398.  
  399. C<< FIN RECHERCHE OBSTRUCTEURS--------------------------------------
  400.  
  401. IF(KIMP.GE.4) THEN
  402. WRITE(6,*) ' FACES K1 K2 ',K1,K2,' NOBS '
  403. $ ,NOBS
  404. ENDIF
  405.  
  406. CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX,NGAX,FF
  407. $ ,KIMP,EXTINC,RAD)
  408. IF(KIMP.GE.4) WRITE(6,*) ' FF ',K1,K2,FF/S1
  409.  
  410. ENDIF
  411.  
  412. ENDIF
  413. C<< FIN CALCUL
  414. C ---------------------------------------------------------------
  415.  
  416. C WRITE(6,*) ' K1 K2 PLIG ',K1,K2,PLIG
  417. PLIG.FACT(K2) = FF/S1
  418.  
  419. ENDIF
  420. C... fin face K2 ..................................................
  421.  
  422.  
  423. C... face inverse de K2 (cas des coques) ..........................
  424.  
  425. IF (ICOQ.AND.KCOQ(JTYP)) THEN
  426.  
  427. K2=K2 + 1
  428.  
  429. KVU=1
  430. FF=0.D0
  431.  
  432. C UTILISATION DE LA RECIPROCITE
  433. C
  434. IF(K1.GT.K2) THEN
  435. S2=PSUR.FACT(K2)
  436. PCOL=LFACT(K2)
  437.  
  438. PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1)
  439.  
  440. ELSE
  441.  
  442. C>> TEST K1=K2 ET VISIBILITE A PRIORI
  443.  
  444. C------------------------------------------------------------------
  445.  
  446. C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0
  447. C
  448. IF(K1.EQ.K2) THEN
  449. IF(KIMP.GE.4) THEN
  450. WRITE(6,*) ' A1 ',A1(1,1),A1(2,1)
  451. WRITE(6,*) ' A1 ',A1(1,2),A1(2,2)
  452. WRITE(6,*) ' U1 ',U1(1),U1(2)
  453. ENDIF
  454. IF (ABS(U1(1)).LT.EPS) KVU=0
  455. IF (U1(1).GE.EPS) KVU=0
  456.  
  457. IF(KVU.NE.0) THEN
  458. DO 514 K = 1,NES
  459. A2(K,1) = A1(K,1)
  460. A2(K,2) = A1(K,2)
  461. 514 CONTINUE
  462. ENDIF
  463.  
  464. C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE
  465. C UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z
  466. C ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z)
  467.  
  468.  
  469. ELSE
  470.  
  471. DO 512 K = 1,NES
  472. U2(K)=-U2(K)
  473. XX1= A2(K,1)
  474. A2(K,1) = A2(K,2) + ECOQ*U2(K)
  475. A2(K,2) = XX1 + ECOQ*U2(K)
  476. 512 CONTINUE
  477.  
  478. C... calcul du symetrique
  479. AA2(1,1)= -A2(1,2)
  480. AA2(2,1)= A2(2,2)
  481. AA2(1,2)= -A2(1,1)
  482. AA2(2,2)= A2(2,1)
  483. CALL KNORM2(NES,AA2,NS,UU2,SS2)
  484.  
  485. C... visibilité a priori
  486. CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU)
  487. CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS)
  488.  
  489. IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS ',K1
  490. $ ,K2,IVU,IVUS
  491. IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0
  492.  
  493. ENDIF
  494. C ---------------------------------------------------------------
  495. C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI
  496. IF(KIMP.GE.4.AND.KVU.NE.0) THEN
  497. WRITE(6,*) ' FACES VISIBLES ',K1,K2
  498. ENDIF
  499.  
  500. C>> CALCUL --------------------------------------------------------
  501. C.. option convexe
  502. IF(KACHE.EQ.0) THEN
  503. IF(KVU.NE.0) THEN
  504.  
  505. CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP,EXTINC
  506. $ ,RAD)
  507. IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE ',K1
  508. $ ,K2,FF/S1
  509.  
  510. ENDIF
  511. C.. option cache
  512. ELSE
  513. IF(KVU.NE.0) THEN
  514.  
  515. C>> les obstructeurs potentiels sont les memes que pour la face K2
  516.  
  517. IF(KIMP.GE.4) THEN
  518. WRITE(6,*) ' FACES K1 K2 ',K1,K2
  519. $ ,' NOBS ',NOBS
  520. ENDIF
  521.  
  522. CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX,NGAX
  523. $ ,FF,KIMP,EXTINC,RAD)
  524.  
  525.  
  526. IF(KIMP.GE.4) WRITE(6,*) ' FF ',K1,K2,FF
  527. $ /S1
  528.  
  529. ENDIF
  530.  
  531. ENDIF
  532. C<< FIN CALCUL ------------------------------------------------------
  533.  
  534. C WRITE(6,*) ' K1 K2 PLIG ',K1,K2,PLIG
  535. PLIG.FACT(K2) = FF/S1
  536.  
  537. ENDIF
  538.  
  539. ENDIF
  540.  
  541. C... fin face inverse de K2 (cas des coques) .........................
  542.  
  543. 210 CONTINUE
  544.  
  545. IF (ICOQ.AND.KCOQ(JTYP)) THEN
  546. NELT2 = NELT2 + 2 * NEL2
  547. ELSE
  548. NELT2 = NELT2 + NEL2
  549. ENDIF
  550.  
  551. 200 CONTINUE
  552.  
  553. C ---FIN BOUCLE FACE 2-----------------------------------------
  554.  
  555.  
  556.  
  557.  
  558. C*** fin traitement de la face K1 ***********************************
  559.  
  560. C*** traitement de la face inverse de K1*****************************
  561. IF (ICOQ.AND.KCOQ(ITYP)) THEN
  562.  
  563. K1=K1+1
  564. PLIG=LFACT(K1)
  565.  
  566. DO 122 K = 1,NES
  567. U1(K) =-U1(K)
  568. XX1=A1(K,1)
  569. A1(K,1) = A1(K,2) + ECOQ * U1(K)
  570. A1(K,2) = XX1 + ECOQ * U1(K)
  571. 122 CONTINUE
  572.  
  573. PSUR.FACT(K1)=S1
  574.  
  575. C>> BOUCLE FACE 2
  576. C -------------------------------------------------------------
  577.  
  578. NELT2= 0
  579. DO 2000 JTYP = 1,NTYP
  580. IMODEL = MYMOD.KMODEL(JTYP)
  581. SEGACT IMODEL
  582. IPT2 = IMAMOD
  583. SEGDES IMODEL
  584.  
  585. NSGEO2 = IPT2.NUM(/1)
  586. NSPA2=1
  587. IF(ICOQ.AND.KQUAD(JTYP)) NSPA2=2
  588. NEL2 = IPT2.NUM(/2)
  589. DO 2100 I2 = 1,NEL2
  590.  
  591. IF (ICOQ.AND.KCOQ(JTYP))THEN
  592. K2 = 2*I2-1 + NELT2
  593. ELSE
  594. K2 = I2 + NELT2
  595. ENDIF
  596.  
  597. C... face K2 .....................................................
  598.  
  599. KVU=1
  600. FF=0.D0
  601.  
  602. C.. UTILISATION DE LA RECIPROCITE
  603. C
  604. IF(K1.GT.K2) THEN
  605. S2=PSUR.FACT(K2)
  606. PCOL=LFACT(K2)
  607.  
  608. PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1)
  609.  
  610. ELSE
  611. C.. TEST K1=K2 ET VISIBILITE A PRIORI
  612.  
  613. C------------------------------------------------------------------
  614.  
  615. C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0
  616. C
  617. IF(K1.EQ.K2) THEN
  618. IF(KIMP.GE.4) THEN
  619. WRITE(6,*) ' A1 ',A1(1,1),A1(2,1)
  620. WRITE(6,*) ' A1 ',A1(1,2),A1(2,2)
  621. WRITE(6,*) ' U1 ',U1(1),U1(2)
  622. ENDIF
  623. IF (ABS(U1(1)).LT.EPS) KVU=0
  624. IF (U1(1).GE.EPS) KVU=0
  625.  
  626. IF(KVU.NE.0) THEN
  627. DO 2140 K = 1,NES
  628. A2(K,1) = A1(K,1)
  629. A2(K,2) = A1(K,2)
  630. 2140 CONTINUE
  631. ENDIF
  632.  
  633. C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE
  634. C UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z
  635. C ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z)
  636.  
  637. ELSE
  638. DO 2110 IS = 1,NSGEO2,NSPA2
  639. LS=IS
  640. IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  641. IREF = (IDIM+1)*(IPT2.NUM(IS,I2)-1)
  642. DO 2120 K = 1,NES
  643. C** A2(K,IS) = XCOOR(IREF+K)/RAD
  644. A2(K,LS) = XCOOR(IREF+K)/RAD
  645. 2120 CONTINUE
  646. 2110 CONTINUE
  647. CALL KNORM2(NES,A2,NS,U2,S2)
  648. S2=XPI*S2*ABS(A2(1,1)+A2(1,2))
  649.  
  650. C.. calcul du symetrique
  651. AA2(1,1)= -A2(1,2)
  652. AA2(2,1)= A2(2,2)
  653. AA2(1,2)= -A2(1,1)
  654. AA2(2,2)= A2(2,1)
  655. CALL KNORM2(NES,AA2,NS,UU2,SS2)
  656.  
  657. C.. visibilite a priori
  658. CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU)
  659. CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS)
  660.  
  661. IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS ',K1
  662. $ ,K2,IVU,IVUS
  663. IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0
  664.  
  665. ENDIF
  666. C ---------------------------------------------------------------
  667. C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI
  668. IF(KIMP.GE.4.AND.KVU.NE.0) THEN
  669. WRITE(6,*) ' FACES VISIBLES ',K1,K2
  670. ENDIF
  671.  
  672. C>> CALCUL
  673. C... option convexe
  674. IF(KACHE.EQ.0) THEN
  675. IF(KVU.NE.0) THEN
  676.  
  677. CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP,EXTINC
  678. $ ,RAD)
  679. IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE ',K1
  680. $ ,K2,FF/S1
  681.  
  682. ENDIF
  683. C.. option cache
  684. ELSE
  685. IF(KVU.NE.0) THEN
  686.  
  687. C>> RECHERCHE DES OBSTRUCTEURS POTENTIELS.--------------------------
  688.  
  689. NOBS=0
  690.  
  691. C>> BOUCLE FACE 3---------------------------------------------------
  692. NELT3= 0
  693. DO 3000 KTYP = 1,NTYP
  694. IMODEL = MYMOD.KMODEL(KTYP)
  695. SEGACT IMODEL
  696. IPT3 = IMAMOD
  697. SEGDES IMODEL
  698.  
  699. NSGEO3 = IPT3.NUM(/1)
  700. NSPA3=1
  701. IF(ICOQ.AND.KQUAD(KTYP)) NSPA3=2
  702. NEL3 = IPT3.NUM(/2)
  703. DO 3100 I3 = 1,NEL3
  704. K3 = I3+ NELT3
  705.  
  706. IF(K3.NE.K1.AND.K3.NE.K2) THEN
  707.  
  708. DO 3110 IS = 1,NSGEO3,NSPA3
  709. LS=IS
  710. IF(ICOQ.AND.KQUAD(ITYP)) LS
  711. $ =(IS+1)/2
  712. IREF = (IDIM+1)*(IPT3.NUM(IS
  713. $ ,I3)-1)
  714. DO 3120 K = 1,NES
  715. C** A3(K,IS) = XCOOR(IREF+K)/RAD
  716. A3(K,LS) = XCOOR(IREF+K)
  717. $ /RAD
  718. 3120 CONTINUE
  719. 3110 CONTINUE
  720.  
  721. ZMIN2=DMIN1(A2(2,1),A2(2,2))
  722. ZMAX2=DMAX1(A2(2,1),A2(2,2))
  723. RMAX2=DMAX1(A2(1,1),A2(1,2))
  724. ZTMIN=DMIN1(ZMIN1,ZMIN2)
  725. ZTMAX=DMAX1(ZMAX1,ZMAX2)
  726. RMAX =DMAX1(RMAX1,RMAX2)
  727.  
  728. IF(DMAX1(A3(2,1),A3(2,2)).LE
  729. $ .ZTMIN) THEN
  730. KOBS=0
  731. ELSEIF(DMIN1(A3(2,1),A3(2,2)).GE
  732. $ .ZTMAX) THEN
  733. KOBS=0
  734. ELSEIF(DMIN1(A3(1,1),A3(1,2)).GE
  735. $ .RMAX) THEN
  736. KOBS=0
  737. ELSE
  738. KOBS=1
  739. OBS(1,2*NOBS+1)=A3(1,1)
  740. OBS(2,2*NOBS+1)=A3(2,1)
  741. OBS(1,2*NOBS+2)=A3(1,2)
  742. OBS(2,2*NOBS+2)=A3(2,2)
  743. NOBS=NOBS+1
  744.  
  745. ENDIF
  746. ENDIF
  747.  
  748. 3100 CONTINUE
  749. C
  750. NELT3=NELT3+NEL3
  751. C
  752. 3000 CONTINUE
  753.  
  754. C<< FIN RECHERCHE OBSTRUCTEURS--------------------------------------
  755.  
  756. IF(KIMP.GE.4) THEN
  757. WRITE(6,*) ' FACES K1 K2 ',K1,K2
  758. $ ,' NOBS ',NOBS
  759. ENDIF
  760.  
  761. CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX,NGAX
  762. $ ,FF,KIMP,EXTINC,RAD)
  763. IF(KIMP.GE.4) WRITE(6,*) ' FF ',K1,K2,FF
  764. $ /S1
  765.  
  766. ENDIF
  767.  
  768. ENDIF
  769. C<< FIN CALCUL
  770. C ---------------------------------------------------------------
  771.  
  772. C WRITE(6,*) ' K1 K2 PLIG ',K1,K2,PLIG
  773. PLIG.FACT(K2) = FF/S1
  774.  
  775. ENDIF
  776. C... fin face K2 ..................................................
  777.  
  778.  
  779. C... face inverse de K2 (cas des coques) ..........................
  780.  
  781. IF (ICOQ.AND.KCOQ(JTYP)) THEN
  782.  
  783. K2=K2 + 1
  784.  
  785. KVU=1
  786. FF=0.D0
  787.  
  788. C UTILISATION DE LA RECIPROCITE
  789. C
  790. IF(K1.GT.K2) THEN
  791. S2=PSUR.FACT(K2)
  792. PCOL=LFACT(K2)
  793.  
  794. PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1)
  795.  
  796. ELSE
  797. C**
  798.  
  799. C>> TEST K1=K2 ET VISIBILITE A PRIORI
  800.  
  801. C------------------------------------------------------------------
  802.  
  803. C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0
  804. C
  805. IF(K1.EQ.K2) THEN
  806. IF(KIMP.GE.4) THEN
  807. WRITE(6,*) ' A1 ',A1(1,1),A1(2,1)
  808. WRITE(6,*) ' A1 ',A1(1,2),A1(2,2)
  809. WRITE(6,*) ' U1 ',U1(1),U1(2)
  810. ENDIF
  811. IF (ABS(U1(1)).LT.EPS) KVU=0
  812. IF (U1(1).GE.EPS) KVU=0
  813.  
  814. IF(KVU.NE.0) THEN
  815. DO 5140 K = 1,NES
  816. A2(K,1) = A1(K,1)
  817. A2(K,2) = A1(K,2)
  818. 5140 CONTINUE
  819. ENDIF
  820.  
  821. C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE
  822. C UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z
  823. C ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z)
  824.  
  825.  
  826. ELSE
  827.  
  828. DO 5120 K = 1,NES
  829. U2(K)=-U2(K)
  830. XX1=A2(K,1)
  831. A2(K,1) = A2(K,2) + ECOQ*U2(K)
  832. A2(K,2) = XX1 + ECOQ*U2(K)
  833. 5120 CONTINUE
  834.  
  835. C... calcul du symetrique
  836. AA2(1,1)= -A2(1,2)
  837. AA2(2,1)= A2(2,2)
  838. AA2(1,2)= -A2(1,1)
  839. AA2(2,2)= A2(2,1)
  840. CALL KNORM2(NES,AA2,NS,UU2,SS2)
  841.  
  842. C... visibilité a priori
  843. CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU)
  844. CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS)
  845.  
  846. IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS '
  847. $ ,K1,K2,IVU,IVUS
  848. IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0
  849.  
  850. ENDIF
  851. C ---------------------------------------------------------------
  852. C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI
  853. IF(KIMP.GE.4.AND.KVU.NE.0) THEN
  854. WRITE(6,*) ' FACES VISIBLES ',K1,K2
  855. ENDIF
  856.  
  857. C>> CALCUL --------------------------------------------------------
  858. C.. option convexe
  859. IF(KACHE.EQ.0) THEN
  860. IF(KVU.NE.0) THEN
  861.  
  862. CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP
  863. $ ,EXTINC,RAD)
  864. IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE '
  865. $ ,K1,K2,FF/S1
  866.  
  867. ENDIF
  868. C.. option cache
  869. ELSE
  870. IF(KVU.NE.0) THEN
  871.  
  872. C>> les obstructeurs potentiels sont les memes que pour la face K2
  873.  
  874. IF(KIMP.GE.4) THEN
  875. WRITE(6,*) ' FACES K1 K2 ',K1,K2
  876. $ ,' NOBS ',NOBS
  877. ENDIF
  878.  
  879. CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX
  880. $ ,NGAX,FF,KIMP,EXTINC,RAD)
  881.  
  882. IF(KIMP.GE.4) WRITE(6,*) ' FF ',K1,K2
  883. $ ,FF/S1
  884.  
  885. ENDIF
  886.  
  887. ENDIF
  888. C<< FIN CALCUL ------------------------------------------------------
  889.  
  890. C WRITE(6,*) ' K1 K2 PLIG ',K1,K2,PLIG
  891. PLIG.FACT(K2) = FF/S1
  892.  
  893. ENDIF
  894.  
  895. ENDIF
  896.  
  897. C... fin face inverse de K2 (cas des coques) .........................
  898.  
  899. 2100 CONTINUE
  900.  
  901. IF (ICOQ.AND.KCOQ(JTYP)) THEN
  902. NELT2 = NELT2 + 2 * NEL2
  903. ELSE
  904. NELT2 = NELT2 + NEL2
  905. ENDIF
  906.  
  907. 2000 CONTINUE
  908.  
  909. C ---FIN BOUCLE FACE 2-----------------------------------------
  910.  
  911. ENDIF
  912. C*** fin traitement de la face inverse de K1*************************
  913.  
  914. 110 CONTINUE
  915.  
  916.  
  917. IF (ICOQ.AND.KCOQ(ITYP)) THEN
  918. NELT1 = NELT1 + 2 * NEL1
  919. ELSE
  920. NELT1 = NELT1 + NEL1
  921. ENDIF
  922.  
  923. 100 CONTINUE
  924.  
  925. C ---FIN BOUCLE FACE 1-----------------------------------------
  926. C
  927.  
  928.  
  929.  
  930. C Désactivation des maillages élémentaires
  931. DO 920 ITYP = 1,NTYP
  932. IMODEL = MYMOD.KMODEL(ITYP)
  933. SEGACT IMODEL
  934. IPT1 = IMAMOD
  935. SEGDES IPT1
  936. SEGDES IMODEL
  937. 920 CONTINUE
  938. SEGDES MYMOD
  939.  
  940. C>> SURFACES DIMENSIONNEES
  941.  
  942. DO 910 LS=1,NFACE
  943. PSUR.FACT(LS)=PSUR.FACT(LS)*RAD*RAD
  944. 910 CONTINUE
  945. C
  946. IF(KACHE.NE.0) SEGSUP SFOBS
  947. C
  948.  
  949. C>>> NORMALISATION,CACUL DES BILANS ET IMPRESSION
  950. C --------------------------------------------
  951.  
  952. IF(EXTINC.GT.1D-3) THEN
  953. INOR=0
  954. ENDIF
  955.  
  956. CALL KFN(IFACFO,INOR,KIMP)
  957.  
  958. C Traduction puis suppression de IFACFO
  959.  
  960. IF (KIMP.GE.3) THEN
  961. CALL PRFACF(IFACFO)
  962. ENDIF
  963.  
  964. LTITR=1
  965. CALL FFMCHA(MYMOD,INFOEL,IFACFO,ICHFAC,LTITR)
  966.  
  967. SEGACT IFACFO
  968. DO 930 NNEL = 1,LFACT(/1)
  969. LFAC = LFACT(NNEL)
  970. SEGSUP LFAC
  971. 930 CONTINUE
  972. SEGSUP IFACFO
  973.  
  974. SEGSUP STRAV
  975.  
  976. RETURN
  977. END
  978. C///////////
  979. C mettre apres KPRIOR
  980. * WRITE (6,*) 'K1 =',K1,' et K2 =',K2
  981. * WRITE (6,*) 'NES =',NES,'NS =',NS
  982. * WRITE (6,*) 'A1(*,1) :',(A1(K,1),K=1,IDIM)
  983. * WRITE (6,*) 'A1(*,2) :',(A1(K,2),K=1,IDIM)
  984. * WRITE (6,*) 'U1 :',(U1(K),K=1,IDIM+1)
  985. * WRITE (6,*) 'A2(*,1) :',(A2(K,1),K=1,IDIM)
  986. * WRITE (6,*) 'A2(*,2) :',(A2(K,2),K=1,IDIM)
  987. * WRITE (6,*) 'U2 :',(U2(K),K=1,IDIM+1)
  988. * WRITE (6,*) 'IVU =',IVU
  989. * WRITE (6,*) 'AA2(*,1) :',(AA2(K,1),K=1,IDIM)
  990. * WRITE (6,*) 'AA2(*,2) :',(AA2(K,2),K=1,IDIM)
  991. * WRITE (6,*) 'UU2 :',(UU2(K),K=1,IDIM+1)
  992. * WRITE (6,*) 'IVUS =',IVUS
  993. C///////////
  994.  
  995.  
  996.  
  997.  
  998.  
  999.  
  1000.  
  1001.  
  1002.  
  1003.  
  1004.  
  1005.  
  1006.  
  1007.  
  1008.  
  1009.  

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