Télécharger ffaxca.eso

Retour à la liste

Numérotation des lignes :

ffaxca
  1. C FFAXCA SOURCE MB234859 21/11/22 21:15:05 11194
  2. SUBROUTINE FFAXCA (ithr,IDONNEES)
  3. C
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8 (A-H,O-Z)
  6. C
  7. LOGICAL ICOQ,IFACINF
  8. C----------------------------------------------------------------------
  9. C Calcul des facteurs de forme en axisymetrique (appele par FACAXI)
  10. C----------------------------------------------------------------------
  11. -INC CCREEL
  12. -INC SMELEME
  13. -INC SMMODEL
  14.  
  15. -INC PPARAM
  16. -INC CCOPTIO
  17. -INC SMCOORD
  18. POINTEUR MYMOD.MMODEL
  19. C
  20. C SEGMENT POUR LA PARALLELISATION
  21. SEGMENT SPARAL
  22. INTEGER NBTHRD
  23. INTEGER IAFAIR(NBEL2)
  24. INTEGER IMYMOD,ISEGEL,KNPAX,KNGAX,KKACHE,IIFACFO
  25. INTEGER KITYP,KNELT1,IWRKTH
  26. REAL*8 XEXTINC,XRAD
  27. LOGICAL BICOQ
  28. ENDSEGMENT
  29. C
  30. C----------------------------------------------------------------------
  31. SEGMENT , INFOEL
  32. LOGICAL KCOQ(N1),KQUAD(N1)
  33. ENDSEGMENT
  34. C----------------------------------------------------------------------
  35. C FACTEURS DE FORME
  36. C NNBEL1 = NOMBRE DE LIGNES + 1
  37. C NBEL2 = NOMBRE DE COLONNES
  38. C LFACT(NNBEL1) POINTE SUR LE TABLEAU DES SURFACES
  39. C
  40. SEGMENT IFACFO
  41. INTEGER LFACT(NNBEL1)
  42. ENDSEGMENT
  43. SEGMENT LFAC
  44. REAL*8 FACT(NBEL2)
  45. ENDSEGMENT
  46. POINTEUR PSUR.LFAC, PLIG.LFAC, PCOL.LFAC
  47. C----------------------------------------------------------------------
  48. C coordonnees des obstructeurs
  49. SEGMENT SFOBS
  50. REAL*8 OBS(2,NTOBS)
  51. ENDSEGMENT
  52. C----------------------------------------------------------------------
  53. SEGMENT STRAV
  54. REAL*8 A1(NA,NA),A2(NA,NA),A3(NA,NA),AA2(NA,NA)
  55. REAL*8 U1(NA1),U2(NA1),UU2(NA1)
  56. ENDSEGMENT
  57. C----------------------------------------------------------------------
  58. SEGMENT SEGTH
  59. INTEGER SSFOBS(NTHRD)
  60. INTEGER SSTRAV(NTHRD)
  61. ENDSEGMENT
  62. C
  63. NS = 2
  64. EPS = 1D-5
  65. KIMP = IIMPI
  66. NES = IDIM
  67. C
  68. SPARAL = IDONNEES
  69. NBTHR = SPARAL.NBTHRD
  70. MYMOD = SPARAL.IMYMOD
  71. INFOEL = SPARAL.ISEGEL
  72. NPAX = SPARAL.KNPAX
  73. NGAX = SPARAL.KNGAX
  74. KACHE = SPARAL.KKACHE
  75. EXTINC = SPARAL.XEXTINC
  76. IFACFO = SPARAL.IIFACFO
  77. ITYP = SPARAL.KITYP
  78. RAD = SPARAL.XRAD
  79. ICOQ = SPARAL.BICOQ
  80. NELT1 = SPARAL.KNELT1
  81. SEGTH = SPARAL.IWRKTH
  82. SFOBS = SEGTH.SSFOBS(ithr)
  83. STRAV = SEGTH.SSTRAV(ithr)
  84. CC WRITE(*,*) 'ith SFOBS STRAV',ithr,sfobs,strav
  85. C
  86. NNBEL1 = LFACT(/1)
  87. PSUR = LFACT(NNBEL1)
  88. NTYP = MYMOD.KMODEL(/1)
  89. C
  90. IMODEL = MYMOD.KMODEL(ITYP)
  91. IPT1 = IMAMOD
  92. NSGEO1 = IPT1.NUM(/1)
  93. NSPA1=1
  94. IF(ICOQ.AND.KQUAD(ITYP)) NSPA1=2
  95. NEL1 = IPT1.NUM(/2)
  96. C
  97. C Decoupage pour le travail d'ecriture en parallele
  98. IRES=MOD(NEL1,NBTHR)
  99. IF (IRES.EQ.0) THEN
  100. ILON = NEL1 / NBTHR
  101. IDEB = (ithr -1)* ILON + 1
  102. ELSE
  103. IF (ithr.LE.IRES) THEN
  104. ILON = (NEL1 / NBTHR) + 1
  105. IDEB = (ithr -1)* ILON + 1
  106. ELSE
  107. ILON = NEL1 / NBTHR
  108. IDEB = (IRES * (ILON+1)) + (ithr-IRES-1)* ILON + 1
  109. ENDIF
  110. ENDIF
  111. IFIN = IDEB + ILON - 1
  112. C
  113. DO 110 I1 = IDEB,IFIN
  114.  
  115. IF (ICOQ.AND.KCOQ(ITYP)) THEN
  116. K1 = (2*I1-1) + NELT1
  117. ELSE
  118. K1 = I1+ NELT1
  119. ENDIF
  120.  
  121. PLIG = LFACT(K1)
  122.  
  123. C*** traitement de la face K1 ***************************************
  124.  
  125. DO 111 IS = 1,NSGEO1,NSPA1
  126. LS=IS
  127. IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  128. IREF = (IDIM+1)*(IPT1.NUM(IS,I1)-1)
  129. IF(KIMP.GE.4) THEN
  130. WRITE (6,*) 'Noeud numero',IPT1.NUM(IS,I1),'ref :',IREF
  131. ENDIF
  132. DO 112 K = 1,NES
  133. A1(K,LS) = XCOOR(IREF+K)/RAD
  134. 112 CONTINUE
  135. 111 CONTINUE
  136.  
  137. CALL KNORM2(NES,A1,NS,U1,S1)
  138. S1=XPI*S1*ABS(A1(1,1)+A1(1,2))
  139. PSUR.FACT(K1)=S1
  140. ZMIN1= MIN(A1(2,1),A1(2,2))
  141. ZMAX1= MAX(A1(2,1),A1(2,2))
  142. RMAX1= MAX(A1(1,1),A1(1,2))
  143. C
  144. C>> BOUCLE FACE 2
  145. C -----------------------------------------------------------
  146. C
  147. IFACINF = .TRUE.
  148. 10 CONTINUE
  149. C
  150. NELT2= 0
  151. DO 200 JTYP = 1,NTYP
  152. C
  153. IMODEL = MYMOD.KMODEL(JTYP)
  154. IPT2 = IMAMOD
  155. NSGEO2 = IPT2.NUM(/1)
  156. NSPA2=1
  157. IF(ICOQ.AND.KQUAD(JTYP)) NSPA2=2
  158. NEL2 = IPT2.NUM(/2)
  159. C
  160. DO 210 I2 = 1,NEL2
  161. IF (ICOQ.AND.KCOQ(JTYP))THEN
  162. K2 = 2*I2-1 + NELT2
  163. ELSE
  164. K2 = I2 + NELT2
  165. ENDIF
  166.  
  167. C... face K2 .....................................................
  168.  
  169. KVU=1
  170. FF=0.D0
  171.  
  172. C.. UTILISATION DE LA RECIPROCITE
  173. C
  174. IF(K1.GT.K2) THEN
  175. IF (NBTHR.NE.1) THEN
  176. IAFAIR(K1) = K2
  177. ELSE
  178. S2=PSUR.FACT(K2)
  179. PCOL=LFACT(K2)
  180. PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1)
  181. ENDIF
  182. ELSE
  183. C.. TEST K1=K2 ET VISIBILITE A PRIORI
  184.  
  185. C------------------------------------------------------------------
  186.  
  187. C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0
  188. C
  189. IF(K1.EQ.K2) THEN
  190. IF(KIMP.GE.4) THEN
  191. WRITE(6,*) ' A1 ',A1(1,1),A1(2,1)
  192. WRITE(6,*) ' A1 ',A1(1,2),A1(2,2)
  193. WRITE(6,*) ' U1 ',U1(1),U1(2)
  194. ENDIF
  195. IF (ABS(U1(1)).LT.EPS) KVU=0
  196. IF (U1(1).GE.EPS) KVU=0
  197.  
  198. IF(KVU.NE.0) THEN
  199. DO 214 K = 1,NES
  200. A2(K,1) = A1(K,1)
  201. A2(K,2) = A1(K,2)
  202. 214 CONTINUE
  203. ENDIF
  204.  
  205. C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE
  206. C UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z
  207. C ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z)
  208.  
  209. ELSE
  210.  
  211. DO 211 IS = 1,NSGEO2,NSPA2
  212. LS=IS
  213. IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  214. IREF = (IDIM+1)*(IPT2.NUM(IS,I2)-1)
  215. DO 212 K = 1,NES
  216. A2(K,LS) = XCOOR(IREF+K)/RAD
  217. 212 CONTINUE
  218. 211 CONTINUE
  219. CALL KNORM2(NES,A2,NS,U2,S2)
  220. S2=XPI*S2*ABS(A2(1,1)+A2(1,2))
  221.  
  222. C.. calcul du symetrique
  223. AA2(1,1)= -A2(1,2)
  224. AA2(2,1)= A2(2,2)
  225. AA2(1,2)= -A2(1,1)
  226. AA2(2,2)= A2(2,1)
  227. CALL KNORM2(NES,AA2,NS,UU2,SS2)
  228.  
  229. C.. visibilite a priori
  230. CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU)
  231. CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS)
  232.  
  233. IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS ',K1,K2,IVU,IVUS
  234. IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0
  235.  
  236. ENDIF
  237. C ---------------------------------------------------------------
  238. C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI
  239. IF(KIMP.GE.4.AND.KVU.NE.0) THEN
  240. WRITE(6,*) ' FACES VISIBLES ',K1,K2
  241. ENDIF
  242.  
  243. C>> CALCUL
  244. C... option convexe
  245. IF(KACHE.EQ.0) THEN
  246. IF(KVU.NE.0) THEN
  247. CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP,EXTINC,RAD)
  248. IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE ',K1,K2,FF/S1
  249.  
  250. ENDIF
  251. ELSE
  252. C.. option cache
  253. IF(KVU.NE.0) THEN
  254.  
  255. C>> RECHERCHE DES OBSTRUCTEURS POTENTIELS.--------------------------
  256.  
  257. NOBS=0
  258.  
  259. ZMIN2= MIN(A2(2,1),A2(2,2))
  260. ZMAX2= MAX(A2(2,1),A2(2,2))
  261. RMAX2= MAX(A2(1,1),A2(1,2))
  262. ZTMIN= MIN(ZMIN1,ZMIN2)
  263. ZTMAX= MAX(ZMAX1,ZMAX2)
  264. RMAX = MAX(RMAX1,RMAX2)
  265.  
  266. C>> BOUCLE FACE 3---------------------------------------------------
  267. NELT3= 0
  268. DO 300 KTYP = 1,NTYP
  269. IMODEL = MYMOD.KMODEL(KTYP)
  270. IPT3 = IMAMOD
  271. NSGEO3 = IPT3.NUM(/1)
  272. NSPA3=1
  273. IF(ICOQ.AND.KQUAD(KTYP)) NSPA3=2
  274. NEL3 = IPT3.NUM(/2)
  275. DO 310 I3 = 1,NEL3
  276. K3 = I3+ NELT3
  277.  
  278. IF(K3.NE.K1.AND.K3.NE.K2) THEN
  279.  
  280. DO 311 IS = 1,NSGEO3,NSPA3
  281. LS=IS
  282. IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  283. IREF = (IDIM+1)*(IPT3.NUM(IS,I3)-1)
  284. DO 312 K = 1,NES
  285. A3(K,LS) = XCOOR(IREF+K)/RAD
  286. 312 CONTINUE
  287. 311 CONTINUE
  288.  
  289. IF( MAX(A3(2,1),A3(2,2)).LE.ZTMIN) THEN
  290. KOBS=0
  291. ELSEIF(MIN(A3(2,1),A3(2,2)).GE.ZTMAX) THEN
  292. KOBS=0
  293. ELSEIF(MIN(A3(1,1),A3(1,2)).GE.RMAX) THEN
  294. KOBS=0
  295. ELSE
  296. KOBS=1
  297. OBS(1,2*NOBS+1)=A3(1,1)
  298. OBS(2,2*NOBS+1)=A3(2,1)
  299. OBS(1,2*NOBS+2)=A3(1,2)
  300. OBS(2,2*NOBS+2)=A3(2,2)
  301. NOBS=NOBS+1
  302. ENDIF
  303. ENDIF
  304. 310 CONTINUE
  305. C
  306. NELT3=NELT3+NEL3
  307. C
  308. 300 CONTINUE
  309.  
  310. C<< FIN RECHERCHE OBSTRUCTEURS--------------------------------------
  311.  
  312. IF(KIMP.GE.4) THEN
  313. WRITE(6,*) ' FACES K1 K2 ',K1,K2,' NOBS ',NOBS
  314. ENDIF
  315.  
  316. CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX,NGAX,FF
  317. $ ,KIMP,EXTINC,RAD)
  318. IF(KIMP.GE.4) WRITE(6,*) ' FF ',K1,K2,FF/S1
  319.  
  320. ENDIF
  321. ENDIF
  322. C<< FIN CALCUL
  323. C ---------------------------------------------------------------
  324.  
  325. C WRITE(6,*) ' K1=',K1,' K2 PLIG ',K2,PLIG
  326. PLIG.FACT(K2) = FF/S1
  327.  
  328. ENDIF
  329. C... fin face K2 ..................................................
  330.  
  331. C... face inverse de K2 (cas des coques) ..........................
  332.  
  333. IF (ICOQ.AND.KCOQ(JTYP)) THEN
  334. K2=K2 + 1
  335. KVU=1
  336. FF=0.D0
  337.  
  338. C UTILISATION DE LA RECIPROCITE
  339. C
  340. IF(K1.GT.K2) THEN
  341. S2=PSUR.FACT(K2)
  342. PCOL=LFACT(K2)
  343. PLIG.FACT(K2)=(S2/S1)*PCOL.FACT(K1)
  344. ELSE
  345.  
  346. C>> TEST K1=K2 ET VISIBILITE A PRIORI
  347.  
  348. C------------------------------------------------------------------
  349.  
  350. C>> K1=K2 : ELIMINATION DES CAS TRIVIAUX OU F(K1,K2)=0
  351. C
  352. IF(K1.EQ.K2) THEN
  353. IF(KIMP.GE.4) THEN
  354. WRITE(6,*) ' A1 ',A1(1,1),A1(2,1)
  355. WRITE(6,*) ' A1 ',A1(1,2),A1(2,2)
  356. WRITE(6,*) ' U1 ',U1(1),U1(2)
  357. ENDIF
  358. IF (ABS(U1(1)).LT.EPS) KVU=0
  359. IF (U1(1).GE.EPS) KVU=0
  360.  
  361. IF(KVU.NE.0) THEN
  362. DO 514 K = 1,NES
  363. A2(K,1) = A1(K,1)
  364. A2(K,2) = A1(K,2)
  365. 514 CONTINUE
  366. ENDIF
  367.  
  368. C>> K1::K2 : ELIMINATION DES CAS DE NON-VISIBLITE
  369. C UTILISATION DE LA VISIBILITE A PRIORI DANS LE PLAN R-Z
  370. C ON TESTE LA FACE K2 ET SON SYMETRIQUE (-R,Z)
  371.  
  372.  
  373. ELSE
  374.  
  375. DO 512 K = 1,NES
  376. U2(K)=-U2(K)
  377. XX1= A2(K,1)
  378. A2(K,1) = A2(K,2) + ECOQ*U2(K)
  379. A2(K,2) = XX1 + ECOQ*U2(K)
  380. 512 CONTINUE
  381.  
  382. C... calcul du symetrique
  383. AA2(1,1)= -A2(1,2)
  384. AA2(2,1)= A2(2,2)
  385. AA2(1,2)= -A2(1,1)
  386. AA2(2,2)= A2(2,1)
  387. CALL KNORM2(NES,AA2,NS,UU2,SS2)
  388.  
  389. C... visibilite a priori
  390. CALL KPRIOR(NES,NS,NS,A1,A2,U1,U2,IVU)
  391. CALL KPRIOR(NES,NS,NS,A1,AA2,U1,UU2,IVUS)
  392.  
  393. IF(KIMP.GE.4)WRITE(6,*) ' K1 K2 IVU IVUS ',K1,K2,IVU,IVUS
  394. IF(IVU.EQ.0.AND.IVUS.EQ.0) KVU=0
  395.  
  396. ENDIF
  397. C ---------------------------------------------------------------
  398. C<< FIN TEST K1=K2 ET VISIBIITE A PRIORI
  399. IF(KIMP.GE.4.AND.KVU.NE.0) THEN
  400. WRITE(6,*) ' FACES VISIBLES ',K1,K2
  401. ENDIF
  402.  
  403. C>> CALCUL --------------------------------------------------------
  404. C.. option convexe
  405. IF(KACHE.EQ.0) THEN
  406. IF(KVU.NE.0) THEN
  407. CALL KAXC(A1,A2,NPAX,NGAX,FF,KIMP,EXTINC,RAD)
  408. IF(KIMP.GE.4) WRITE(6,*) ' FF CONVEXE ',K1,K2,FF/S1
  409. ENDIF
  410. C.. option cache
  411. ELSE
  412. IF(KVU.NE.0) THEN
  413.  
  414. C>> les obstructeurs potentiels sont les memes que pour la face K2
  415.  
  416. IF(KIMP.GE.4) THEN
  417. WRITE(6,*) ' FACES K1 K2 ',K1,K2,' NOBS ',NOBS
  418. ENDIF
  419.  
  420. CALL KAXK(A1,A2,OBS,NOBS,NTOBS,NPAX,NGAX
  421. $ ,FF,KIMP,EXTINC,RAD)
  422. IF(KIMP.GE.4) WRITE(6,*) ' FF ',K1,K2,FF/S1
  423. ENDIF
  424. ENDIF
  425. C<< FIN CALCUL ------------------------------------------------------
  426.  
  427. C WRITE(6,*) ' K1 K2 PLIG ',K1,K2,PLIG
  428. PLIG.FACT(K2) = FF/S1
  429.  
  430. ENDIF
  431. ENDIF
  432.  
  433. C... fin face inverse de K2 (cas des coques) .........................
  434.  
  435. 210 CONTINUE
  436.  
  437. IF (ICOQ.AND.KCOQ(JTYP)) THEN
  438. NELT2 = NELT2 + 2 * NEL2
  439. ELSE
  440. NELT2 = NELT2 + NEL2
  441. ENDIF
  442.  
  443. 200 CONTINUE
  444.  
  445. C*** traitement de la face inverse de K1*****************************
  446. IF (ICOQ.AND.KCOQ(ITYP).AND.IFACINF) THEN
  447. K1=K1+1
  448. PLIG = LFACT(K1)
  449.  
  450. DO 122 K = 1,NES
  451. U1(K) =-U1(K)
  452. XX1=A1(K,1)
  453. A1(K,1) = A1(K,2) + ECOQ * U1(K)
  454. A1(K,2) = XX1 + ECOQ * U1(K)
  455. 122 CONTINUE
  456.  
  457. PSUR.FACT(K1)=S1
  458. IFACINF = .FALSE.
  459. GOTO 10
  460.  
  461. ENDIF
  462. C*** fin traitement de la face inverse de K1*************************
  463.  
  464. 110 CONTINUE
  465. C
  466. RETURN
  467. END
  468. C///////////
  469. C mettre apres KPRIOR
  470. * WRITE (6,*) 'K1 =',K1,' et K2 =',K2
  471. * WRITE (6,*) 'NES =',NES,'NS =',NS
  472. * WRITE (6,*) 'A1(*,1) :',(A1(K,1),K=1,IDIM)
  473. * WRITE (6,*) 'A1(*,2) :',(A1(K,2),K=1,IDIM)
  474. * WRITE (6,*) 'U1 :',(U1(K),K=1,IDIM+1)
  475. * WRITE (6,*) 'A2(*,1) :',(A2(K,1),K=1,IDIM)
  476. * WRITE (6,*) 'A2(*,2) :',(A2(K,2),K=1,IDIM)
  477. * WRITE (6,*) 'U2 :',(U2(K),K=1,IDIM+1)
  478. * WRITE (6,*) 'IVU =',IVU
  479. * WRITE (6,*) 'AA2(*,1) :',(AA2(K,1),K=1,IDIM)
  480. * WRITE (6,*) 'AA2(*,2) :',(AA2(K,2),K=1,IDIM)
  481. * WRITE (6,*) 'UU2 :',(UU2(K),K=1,IDIM+1)
  482. * WRITE (6,*) 'IVUS =',IVUS
  483. C///////////
  484.  
  485.  

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