Télécharger kalpre.eso

Retour à la liste

Numérotation des lignes :

kalpre
  1. C KALPRE SOURCE CB215821 17/11/30 21:16:36 9639
  2. SUBROUTINE KALPRE (MYMOD,INFOEL,SKFACE,XDEC,NELD)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. C----------------------------------------------------------------------
  6. C Facteurs de forme OPTION CACHE 3D
  7. C Calcul du decoupage des faces et initialisations
  8. C SP appelle par facgen
  9. C entree:
  10. C MYMOD : pointeur de l'objet modèle
  11. C INFOEL : informations concernant le type des éléments des maillages
  12. C SKFACE : pointeur sur l'objet 'faces' pour le calcul des FF
  13. C SHC3D : pointeur sur la surface de projection
  14. C XDEC : paramŠtre pour le decoupage des faces
  15. C NELD : nombre d'éléments virtuels
  16. C (un élément COQ -> 2 facteurs de forme <=> 2 éléments virtuels)
  17. C sortie:
  18. C SKFACE : pointeur sur l'objet 'faces' pour le calcul des FF
  19. C SHC3D : pointeur sur la surface de projection
  20. C----------------------------------------------------------------------
  21. C
  22.  
  23. -INC PPARAM
  24. -INC CCOPTIO
  25. -INC SMCOORD
  26. -INC SMMODEL
  27. -INC SMELEME
  28. -INC TFFOR3D
  29. C
  30. C -----------------------------------
  31. SEGMENT SDECOU
  32. INTEGER KDECOU(MFACE)
  33. ENDSEGMENT
  34. POINTEUR MYMOD.MMODEL
  35. C -----------------------------------
  36. C Stockage d'informations concernant le type des éléments des maillages
  37. SEGMENT ,INFOEL
  38. LOGICAL KCOQ(N1),KQUAD(N1)
  39. ENDSEGMENT
  40. C -----------------------------------
  41. DIMENSION A(3,4),UU(4),IN(3)
  42. DIMENSION A1(3,3),A2(3,3),U1(4),U2(4),G1(3),G2(3)
  43. DIMENSION GG(3,400),SS(400)
  44. LOGICAL ICOQ
  45. C
  46. C--------------------------------------------------------------------
  47. C
  48. C
  49. DMIN=1.E-5
  50. DMAX=1./DMIN
  51. C... nombre max de triangles de subdivision
  52. NPM=20
  53. C... quadratiques
  54. NSPA=1
  55.  
  56. IF(INFOEL.EQ.0) THEN
  57. ICOQ=.FALSE.
  58. ELSE
  59. ICOQ=.TRUE.
  60. SEGACT INFOEL
  61. ENDIF
  62. C
  63. C>>> CREATION DE L'OBJET FACE
  64. C
  65. NFACE = 0
  66. NELD = 0
  67. C
  68. SEGACT MYMOD
  69. MTYP = MYMOD.KMODEL(/1)
  70. SEGACT INFOEL
  71. C
  72. C On va compter le nombre d'éléments virtuels (égal au nombre de
  73. C facteurs de forme par ligne) : NELD , et le nombre de faces : NFACE .
  74. DO 10 ITYP=1,MTYP
  75. IMODE1 = MYMOD.KMODEL(ITYP)
  76. SEGACT IMODE1
  77. IPT1 = IMODE1.IMAMOD
  78. SEGDES IMODE1
  79. SEGACT IPT1
  80. NEL = IPT1.NUM(/2)
  81. NPGEO = IPT1.NUM(/1)
  82. SEGDES IPT1
  83. C
  84. C un élément COQ -> 2 facteurs de formes
  85. IF (KCOQ(ITYP)) NEL=2*NEL
  86. NELD = NELD + NEL
  87. C
  88. C un élément carré -> 2 éléments triangles
  89. IF(NPGEO.EQ.4.OR.NPGEO.EQ.8) NEL = 2*NEL
  90. NFACE=NFACE+NEL
  91. 10 CONTINUE
  92. C
  93. IF (KIMP.GE.1) THEN
  94. WRITE( 6,*) ' NOMBRE TOTAL D ELEMENTS ',NELD,' DE FACES ',NFACE
  95. ENDIF
  96. C
  97. C
  98. MFACE = NFACE
  99. MS = 3
  100. SEGINI SKFACE
  101. SEGINI SDECOU
  102. C
  103. C KEL : INDICE ELEMENT DECOUPAGE (-> numéro global de la facette)
  104. C IEL : INDICE ELEMENT D'ORIGINE (-> numero global de l'élément virtuel)
  105. C Ces indices sont différents car on découpe les quadrangles en triangles
  106. C
  107. KEL = 1
  108. IEL = 0
  109. C
  110. DO 100 ITYP = 1, MTYP
  111. IMODE1 = MYMOD.KMODEL(ITYP)
  112. SEGACT IMODE1
  113. IPT1 = IMODE1.IMAMOD
  114. SEGDES IMODE1
  115. SEGACT IPT1
  116.  
  117. NSGEO = IPT1.NUM(/1)
  118. IF(ICOQ.AND.KQUAD(ITYP)) THEN
  119. NS = NSGEO/2
  120. ELSE
  121. NS=NSGEO
  122. ENDIF
  123. IF(ICOQ.AND.KQUAD(ITYP)) NSPA=2
  124. NEL = IPT1.NUM(/2)
  125. C
  126. C... TRI3 ou TRI6
  127. IF (NS.EQ.3) THEN
  128. *////////
  129. * WRITE (6,*) 'Le maillage ',IPT1,
  130. * # 'est constitué d éléments triangles'
  131. *////////
  132. DO 110 I = 1,NEL
  133. IF (KCOQ(ITYP)) THEN
  134. IEL1 = IEL + (2*I-1)
  135. ELSE
  136. IEL1 = IEL + I
  137. ENDIF
  138. KCORR(KEL) = IEL1
  139.  
  140. DO 111 IS = 1,NSGEO,NSPA
  141. LS=IS
  142. IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  143. C** NUK(IS,KEL) = IPT1.NUM(IS,I)
  144. NUK(LS,KEL) = IPT1.NUM(IS,I)
  145. IREF = (IDIM+1)*(NUK(LS,KEL)-1)
  146. DO 112 K = 1,KES
  147. C** A(K,LS) = XCOOR(IREF+K)
  148. A(K,LS) = XCOOR(IREF+K)
  149. 112 CONTINUE
  150. *////////
  151. * WRITE (6,*) 'noeud numero',IPT1.NUM(IS,I)
  152. * WRITE (6,*) 'Coordonnées :',(A(K,IS),K=1,KES)
  153. *////////
  154. 111 CONTINUE
  155. IF (KIMP.GE.3) THEN
  156. WRITE (6,*) 'La facette ',KEL,'repose sur les points',
  157. # (IPT1.NUM(IS,I),IS=1,NS)
  158. ENDIF
  159. CALL KNORM(KES,A,NS,UU,S(KEL))
  160. DO 113 L = 1,4
  161. U(L,KEL) = UU(L)
  162. 113 CONTINUE
  163. C WRITE(6,*) ' KEL S U ',KEL,S(KEL),(U(K,NEL),K=1,4)
  164.  
  165. KEL = KEL + 1
  166. *////////
  167. IF (KCOQ(ITYP)) THEN
  168. C On remplit KCOOR , NUK , U et S pour l'élément inverse
  169. KCORR(KEL) = IEL1 + 1
  170. NUK(1,KEL) = NUK(3,KEL-1)
  171. NUK(2,KEL) = NUK(2,KEL-1)
  172. NUK(3,KEL) = NUK(1,KEL-1)
  173. DO 114 L = 1,4
  174. U(L,KEL) = - UU(L)
  175. 114 CONTINUE
  176. S(KEL) = S(KEL-1)
  177. KEL = KEL + 1
  178. ENDIF
  179. C
  180. 110 CONTINUE
  181. C
  182. IF (KCOQ(ITYP)) THEN
  183. IEL = IEL + NEL*2
  184. ELSE
  185. IEL = IEL + NEL
  186. ENDIF
  187. C
  188. C... QUA4 ou QUA8
  189. ELSE
  190. *////////
  191. * WRITE (6,*) 'Le maillage ',IPT1,
  192. * # 'est constitué d éléments carrés'
  193. *////////
  194. C
  195. DO 120 I = 1,NEL
  196. C
  197. IF (KCOQ(ITYP)) THEN
  198. IEL1 = IEL + (2*I-1)
  199. ELSE
  200. IEL1 = IEL + I
  201. ENDIF
  202. DO 121 IS = 1,NSGEO,NSPA
  203. LS=IS
  204. IF(ICOQ.AND.KQUAD(ITYP)) LS=(IS+1)/2
  205. IREF = (IDIM+1)*(IPT1.NUM(IS,I)-1)
  206. *////////
  207. * WRITE (6,*) 'noeud numero',IPT1.NUM(IS,I)
  208. * WRITE (6,*) 'IREF =',IREF
  209. *////////
  210. DO 122 K = 1,KES
  211. C** A(K,IS) = XCOOR(IREF+K)
  212. A(K,LS) = XCOOR(IREF+K)
  213. 122 CONTINUE
  214. *////////
  215. * WRITE (6,*) 'noeud numero',IPT1.NUM(IS,I)
  216. * WRITE (6,*) 'Coordonnées :',(A(K,IS),K=1,KES)
  217. *////////
  218. 121 CONTINUE
  219. C
  220. C On regarde comment on va découper le quadrangle
  221. DIA1 = 0.
  222. DO 123 K = 1,KES
  223. DIA1 = DIA1 + (A(K,1)-A(K,3))**2
  224. 123 CONTINUE
  225. DIA2 = 0.
  226. DO 124 K = 1,KES
  227. DIA2 = DIA2 + (A(K,2)-A(K,4))**2
  228. 124 CONTINUE
  229. *//////
  230. * WRITE(6,*) ' DIA1 DIA2 ',DIA1,DIA2
  231. *//////
  232.  
  233. C stockage des connectivités des traiangles
  234.  
  235. IF(ICOQ.AND.KQUAD(ITYP)) THEN
  236.  
  237. * cas des qua8
  238.  
  239. IF (DIA1.LE.DIA2) THEN
  240. NUK(1,KEL) = IPT1.NUM(1,I)
  241. NUK(2,KEL) = IPT1.NUM(3,I)
  242. NUK(3,KEL) = IPT1.NUM(5,I)
  243. KEL = KEL + 1
  244. IF (KCOQ(ITYP)) THEN
  245. NUK(1,KEL) = IPT1.NUM(5,I)
  246. NUK(2,KEL) = IPT1.NUM(3,I)
  247. NUK(3,KEL) = IPT1.NUM(1,I)
  248. KEL = KEL + 1
  249. ENDIF
  250. NUK(1,KEL) = IPT1.NUM(5,I)
  251. NUK(2,KEL) = IPT1.NUM(7,I)
  252. NUK(3,KEL) = IPT1.NUM(1,I)
  253. KEL = KEL + 1
  254. IF (KCOQ(ITYP)) THEN
  255. NUK(1,KEL) = IPT1.NUM(1,I)
  256. NUK(2,KEL) = IPT1.NUM(7,I)
  257. NUK(3,KEL) = IPT1.NUM(5,I)
  258. KEL = KEL + 1
  259. ENDIF
  260. ELSE
  261. NUK(1,KEL) = IPT1.NUM(3,I)
  262. NUK(2,KEL) = IPT1.NUM(5,I)
  263. NUK(3,KEL) = IPT1.NUM(7,I)
  264. KEL = KEL + 1
  265. IF (KCOQ(ITYP)) THEN
  266. NUK(1,KEL) = IPT1.NUM(7,I)
  267. NUK(2,KEL) = IPT1.NUM(5,I)
  268. NUK(3,KEL) = IPT1.NUM(3,I)
  269. KEL = KEL + 1
  270. ENDIF
  271. NUK(1,KEL) = IPT1.NUM(7,I)
  272. NUK(2,KEL) = IPT1.NUM(1,I)
  273. NUK(3,KEL) = IPT1.NUM(3,I)
  274. KEL = KEL + 1
  275. IF (KCOQ(ITYP)) THEN
  276. NUK(1,KEL) = IPT1.NUM(3,I)
  277. NUK(2,KEL) = IPT1.NUM(1,I)
  278. NUK(3,KEL) = IPT1.NUM(7,I)
  279. KEL = KEL + 1
  280. ENDIF
  281. ENDIF
  282.  
  283. ELSE
  284.  
  285. * cas des qua4
  286.  
  287. IF (DIA1.LE.DIA2) THEN
  288. NUK(1,KEL) = IPT1.NUM(1,I)
  289. NUK(2,KEL) = IPT1.NUM(2,I)
  290. NUK(3,KEL) = IPT1.NUM(3,I)
  291. KEL = KEL + 1
  292. IF (KCOQ(ITYP)) THEN
  293. NUK(1,KEL) = IPT1.NUM(3,I)
  294. NUK(2,KEL) = IPT1.NUM(2,I)
  295. NUK(3,KEL) = IPT1.NUM(1,I)
  296. KEL = KEL + 1
  297. ENDIF
  298. NUK(1,KEL) = IPT1.NUM(3,I)
  299. NUK(2,KEL) = IPT1.NUM(4,I)
  300. NUK(3,KEL) = IPT1.NUM(1,I)
  301. KEL = KEL + 1
  302. IF (KCOQ(ITYP)) THEN
  303. NUK(1,KEL) = IPT1.NUM(1,I)
  304. NUK(2,KEL) = IPT1.NUM(4,I)
  305. NUK(3,KEL) = IPT1.NUM(3,I)
  306. KEL = KEL + 1
  307. ENDIF
  308. ELSE
  309. NUK(1,KEL) = IPT1.NUM(2,I)
  310. NUK(2,KEL) = IPT1.NUM(3,I)
  311. NUK(3,KEL) = IPT1.NUM(4,I)
  312. KEL = KEL + 1
  313. IF (KCOQ(ITYP)) THEN
  314. NUK(1,KEL) = IPT1.NUM(4,I)
  315. NUK(2,KEL) = IPT1.NUM(3,I)
  316. NUK(3,KEL) = IPT1.NUM(2,I)
  317. KEL = KEL + 1
  318. ENDIF
  319. NUK(1,KEL) = IPT1.NUM(4,I)
  320. NUK(2,KEL) = IPT1.NUM(1,I)
  321. NUK(3,KEL) = IPT1.NUM(2,I)
  322. KEL = KEL + 1
  323. IF (KCOQ(ITYP)) THEN
  324. NUK(1,KEL) = IPT1.NUM(2,I)
  325. NUK(2,KEL) = IPT1.NUM(1,I)
  326. NUK(3,KEL) = IPT1.NUM(4,I)
  327. KEL = KEL + 1
  328. ENDIF
  329. ENDIF
  330.  
  331. ENDIF
  332.  
  333. C
  334. DO 125 IL = 1,2
  335. C
  336. IF (KCOQ(ITYP)) THEN
  337. INDICE = KEL - IL*2
  338. ELSE
  339. INDICE = KEL - IL
  340. ENDIF
  341. KCORR(INDICE) = IEL1
  342. C
  343. DO 126 IS = 1,3
  344. IREF = (IDIM+1)*(NUK(IS,INDICE)-1)
  345. DO 127 K = 1,KES
  346. A(K,IS) = XCOOR(IREF+K)
  347. C WRITE(6,*) 'A ',K,IS,A(K,IS)
  348. 127 CONTINUE
  349. 126 CONTINUE
  350. CALL KNORM(KES,A,3,UU,S(INDICE))
  351. DO 128 L = 1,4
  352. U(L,INDICE) = UU(L)
  353. 128 CONTINUE
  354. CC WRITE(6,*) ' KEL S U ',INDICE,S(INDICE),(U(K,INDICE),K=1,4)
  355.  
  356. 125 CONTINUE
  357. C
  358. IF (KCOQ(ITYP)) THEN
  359. C On remplit KCOOR , U et S pour les éléments inverses
  360. DO 135 IL = 1,2
  361. INDICE = KEL - IL*2 +1
  362. KCORR(INDICE) = IEL1 +1
  363. DO 138 L = 1,4
  364. U(L,INDICE) = - U(L,INDICE-1)
  365. 138 CONTINUE
  366. S(INDICE)=S(INDICE-1)
  367. 135 CONTINUE
  368. ENDIF
  369. C
  370. 120 CONTINUE
  371. C
  372. IF (KCOQ(ITYP)) THEN
  373. IEL = IEL + NEL*2
  374. ELSE
  375. IEL = IEL + NEL
  376. ENDIF
  377. C
  378. C
  379. ENDIF
  380. C
  381. SEGDES IPT1
  382. C
  383. 100 CONTINUE
  384. C
  385. C On a pris toute l'information nécessaire concernant le maillage
  386. SEGDES MYMOD
  387. SEGDES INFOEL
  388. *////////
  389. IF(KIMP.GE.2) THEN
  390. WRITE (6,*) 'Subdivision'
  391. ENDIF
  392. *////////
  393. C-----------------------------------------------------------
  394.  
  395. DO 200 K1=1,MFACE
  396.  
  397. DO 202 K = 1,KES
  398. G1(K)=0.
  399. DO 201 ISS = 1,MS
  400. IREF = (IDIM+1)*(NUK(ISS,K1)-1)
  401. A1(K,ISS) = XCOOR(IREF+K)
  402. G1(K) = G1(K) + A1(K,ISS)
  403. 201 CONTINUE
  404. G1(K)=G1(K)/3.
  405. 202 CONTINUE
  406. DO 203 K=1,KES+1
  407. U1(K) = U(K,K1)
  408. 203 CONTINUE
  409. C
  410. IF (XDEC.GE.0.01) THEN
  411. DK1 = DMAX
  412. DO 400 K2 = 1,MFACE
  413. C
  414. DO 412 K=1,KES+1
  415. U2(K) = U(K,K2)
  416. 412 CONTINUE
  417.  
  418. DO 212 K = 1,KES
  419. G2(K)=0.
  420. DO 211 ISS = 1,MS
  421. IREF = (IDIM+1)*(NUK(ISS,K2)-1)
  422. A2(K,ISS) = XCOOR(IREF+K)
  423. G2(K)=G2(K)+A2(K,ISS)
  424. 211 CONTINUE
  425. G2(K)=G2(K)/3.
  426. 212 CONTINUE
  427.  
  428. C>>> VISIBILITE A PRIORI
  429. C
  430.  
  431. CALL KPRIOR(KES,MS,MS,A1,A2,U1,U2,KVU)
  432. IF (KVU.NE.0) THEN
  433. D1=G1(1)-G2(1)
  434. D2=G1(2)-G2(2)
  435. D3=G1(3)-G2(3)
  436. DKK1 = SQRT(D1*D1+D2*D2+D3*D3)
  437. C WRITE(6,*) ' DKK1 CK1 ',DKK1,CK1
  438. IF(DKK1.GE.DMIN) THEN
  439. CK1 = ABS(U1(1)*D1+ U1(2)*D2 + U1(3)*D3)/DKK1
  440. IF(CK1.GE.0.01) THEN
  441. DK1 =MIN(DK1,DKK1)
  442. ENDIF
  443. ENDIF
  444. ENDIF
  445. 400 CONTINUE
  446. DR = DK1/XDEC
  447. ELSE
  448. DR = DMAX
  449. ENDIF
  450.  
  451. C
  452. 901 CONTINUE
  453. CALL KCREPA(DR,A1,S(K1),KES,MS,NPAT,GG,SS)
  454.  
  455. IF(NPAT.GT.NPM) THEN
  456. DR=DR*2.
  457.  
  458. IF (KIMP.GE.2) THEN
  459. WRITE(6,900) K1,NPAT
  460. 900 FORMAT(1X,' FACE ',I4,' : SUBDIVISION EXCESSIVE ',I4)
  461. ENDIF
  462. GOTO 901
  463. ENDIF
  464. NPATCH = NPAT
  465. KDECOU(K1)=NPAT
  466. MSP = KES
  467. SEGINI IPATCH
  468.  
  469. DO 130 IP = 1,NPAT
  470. SP(IP) = SS(IP)/S(K1)
  471. DO 129 K=1,KES
  472. GP(K,IP) = GG(K,IP)
  473. 129 CONTINUE
  474. C WRITE(6,*) ' SP GP ',SP(IP),GP(1,IP),GP(2,IP),GP(3,IP)
  475. 130 CONTINUE
  476.  
  477. KPATCH(K1) = IPATCH
  478. SEGDES IPATCH
  479.  
  480. 200 CONTINUE
  481. C-----------------------------------------------------------
  482.  
  483. IF(KIMP.GE.2) THEN
  484. WRITE(6,*)
  485. WRITE(6,*) 'NOMBRE DE SUBDIVISIONS PAR FACE '
  486. WRITE(6,1000) (KDECOU(I),I=1,MFACE)
  487. 1000 FORMAT(1X,10(I4))
  488. ENDIF
  489. C
  490. SEGDES SKFACE
  491. SEGSUP SDECOU
  492. C
  493. RETURN
  494. END
  495.  
  496.  
  497.  
  498.  
  499.  
  500.  
  501.  
  502.  
  503.  
  504.  
  505.  
  506.  

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