Télécharger kdom2.eso

Retour à la liste

Numérotation des lignes :

  1. C KDOM2 SOURCE KK2000 14/04/10 21:15:11 8032
  2. SUBROUTINE KDOM2(MELEMQ)
  3. C
  4. C************************************************************************
  5. C
  6. C PROJET : CASTEM 2000
  7. C
  8. C NOM : KDOM2
  9. C
  10. C DESCRIPTION : Subroutine called by KDOM1
  11. C In the case of EULER model, we change the
  12. C position of the QUAF points
  13. C
  14. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI)
  15. C
  16. C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF
  17. C
  18. C************************************************************************
  19. C
  20. C
  21.  
  22. IMPLICIT INTEGER(I-N)
  23. IMPLICIT REAL*8(A-H,O-Z)
  24.  
  25. INTEGER MELEMQ, NBSOUS, ISOUS, NBEL, IEL, NN1, NN2, NN3, I1
  26. & ,NN4,NN5,NN6,NN7,NN8,NN9,NN10,NN11,NN12,NN13,NN14,NN15
  27. & ,NN16,NN17,NN18,NN19,NN20,NN21,NN22,NN23,NN24,NN25,NN26
  28. & ,NN27
  29. REAL*8 XCEN(3),VOL,X1(3),X2(3),X3(3),X4(3),X5(3),X6(3)
  30. & ,VOL1,VOL2,VOL3,VOL4,VOL5,VOL6
  31. LOGICAL LOSEG3, LOTRI7, LOQUA9, LOTE15, LOPY19, LOPR21, LOCU27
  32. C
  33. -INC CCOPTIO
  34. -INC SMELEME
  35. -INC SMCOORD
  36. C
  37. C Elements allowed are SEG3, TRI7, QUA9, TE15, 'PY19', 'PR21', 'CU27'
  38. C ITYPEL 3 7 11 35 36 34 33
  39. C They can be referred just once
  40. C
  41. LOSEG3=.FALSE.
  42. LOTRI7=.FALSE.
  43. LOQUA9=.FALSE.
  44. LOTE15=.FALSE.
  45. LOPY19=.FALSE.
  46. LOPR21=.FALSE.
  47. LOCU27=.FALSE.
  48. C
  49. MELEME=MELEMQ
  50. SEGACT MELEME
  51. NBSOUS=MELEME.LISOUS(/1)
  52. IF(NBSOUS .EQ. 0) NBSOUS=1
  53. DO ISOUS=1,NBSOUS,1
  54. IF(NBSOUS .NE. 1)THEN
  55. IPT1=MELEME.LISOUS(ISOUS)
  56. SEGACT IPT1
  57. ELSE
  58. IPT1=MELEME
  59. ENDIF
  60. NBEL=IPT1.NUM(/2)
  61. C
  62. IF(IPT1.ITYPEL .EQ. 3)THEN
  63. IF(LOSEG3)THEN
  64. C Elt already referred
  65. WRITE(IOIMP,*) 'Subroutine kdom2'
  66. WRITE(IOIMP,*) 'Mesh type not recognized'
  67. CALL ERREUR(5)
  68. ENDIF
  69. LOSEG3=.TRUE.
  70. C
  71. C********** SEG3
  72. C
  73. DO IEL=1,NBEL,1
  74. NN1=IPT1.NUM(1,IEL)
  75. NN2=IPT1.NUM(2,IEL)
  76. NN3=IPT1.NUM(3,IEL)
  77. DO I1=1,IDIM
  78. XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0*
  79. & (XCOOR((NN3-1)*(IDIM+1)+I1)+
  80. & XCOOR((NN1-1)*(IDIM+1)+I1))
  81. ENDDO
  82. ENDDO
  83. ELSEIF(IPT1.ITYPEL .EQ. 7)THEN
  84. C
  85. C********** TRI7
  86. C
  87. IF(LOTRI7)THEN
  88. C Elt already referred
  89. WRITE(IOIMP,*) 'Subroutine kdom2'
  90. WRITE(IOIMP,*) 'Mesh type not recognized'
  91. CALL ERREUR(5)
  92. ENDIF
  93. LOTRI7=.TRUE.
  94. C
  95. DO IEL=1,NBEL,1
  96. NN1=IPT1.NUM(1,IEL)
  97. NN2=IPT1.NUM(2,IEL)
  98. NN3=IPT1.NUM(3,IEL)
  99. NN4=IPT1.NUM(4,IEL)
  100. NN5=IPT1.NUM(5,IEL)
  101. NN6=IPT1.NUM(6,IEL)
  102. NN7=IPT1.NUM(7,IEL)
  103. DO I1=1,IDIM
  104. XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0*
  105. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  106. & XCOOR((NN3-1)*(IDIM+1)+I1))
  107. XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0*
  108. & (XCOOR((NN3-1)*(IDIM+1)+I1)+
  109. & XCOOR((NN5-1)*(IDIM+1)+I1))
  110. XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0*
  111. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  112. & XCOOR((NN5-1)*(IDIM+1)+I1))
  113. XCOOR((NN7-1)*(IDIM+1)+I1)=
  114. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  115. & XCOOR((NN3-1)*(IDIM+1)+I1)+
  116. & XCOOR((NN5-1)*(IDIM+1)+I1))/3.0D0
  117. ENDDO
  118. ENDDO
  119. ELSEIF(IPT1.ITYPEL .EQ. 11)THEN
  120. C
  121. C********** QUA9
  122. C
  123. IF(LOQUA9)THEN
  124. C Elt already referred
  125. WRITE(IOIMP,*) 'Subroutine kdom2'
  126. WRITE(IOIMP,*) 'Mesh type not recognized'
  127. CALL ERREUR(5)
  128. ENDIF
  129. LOQUA9=.TRUE.
  130. C
  131. DO IEL=1,NBEL,1
  132. NN1=IPT1.NUM(1,IEL)
  133. NN2=IPT1.NUM(2,IEL)
  134. NN3=IPT1.NUM(3,IEL)
  135. NN4=IPT1.NUM(4,IEL)
  136. NN5=IPT1.NUM(5,IEL)
  137. NN6=IPT1.NUM(6,IEL)
  138. NN7=IPT1.NUM(7,IEL)
  139. NN8=IPT1.NUM(8,IEL)
  140. NN9=IPT1.NUM(9,IEL)
  141. C
  142. C**************** Centre de l'elt
  143. C
  144. CALL KDOM7(NN1,NN3,NN5,NN7,XCEN)
  145. IF(IERR.NE.0)GOTO 9999
  146. DO I1=1,IDIM
  147. XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0*
  148. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  149. & XCOOR((NN3-1)*(IDIM+1)+I1))
  150. XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0*
  151. & (XCOOR((NN3-1)*(IDIM+1)+I1)+
  152. & XCOOR((NN5-1)*(IDIM+1)+I1))
  153. XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0*
  154. & (XCOOR((NN7-1)*(IDIM+1)+I1)+
  155. & XCOOR((NN5-1)*(IDIM+1)+I1))
  156. XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0*
  157. & (XCOOR((NN7-1)*(IDIM+1)+I1)+
  158. & XCOOR((NN1-1)*(IDIM+1)+I1))
  159. XCOOR((NN9-1)*(IDIM+1)+I1)=XCEN(I1)
  160. ENDDO
  161. ENDDO
  162. ELSEIF(IPT1.ITYPEL .EQ. 35)THEN
  163. C
  164. C********** TE15
  165. C
  166. IF(LOTE15)THEN
  167. C Elt already referred
  168. WRITE(IOIMP,*) 'Subroutine kdom2'
  169. WRITE(IOIMP,*) 'Mesh type not recognized'
  170. CALL ERREUR(5)
  171. ENDIF
  172. LOTE15=.TRUE.
  173. C
  174. DO IEL=1,NBEL,1
  175. NN1=IPT1.NUM(1,IEL)
  176. NN2=IPT1.NUM(2,IEL)
  177. NN3=IPT1.NUM(3,IEL)
  178. NN4=IPT1.NUM(4,IEL)
  179. NN5=IPT1.NUM(5,IEL)
  180. NN6=IPT1.NUM(6,IEL)
  181. NN7=IPT1.NUM(7,IEL)
  182. NN8=IPT1.NUM(8,IEL)
  183. NN9=IPT1.NUM(9,IEL)
  184. NN10=IPT1.NUM(10,IEL)
  185. NN11=IPT1.NUM(11,IEL)
  186. NN12=IPT1.NUM(12,IEL)
  187. NN13=IPT1.NUM(13,IEL)
  188. NN14=IPT1.NUM(14,IEL)
  189. NN15=IPT1.NUM(15,IEL)
  190. DO I1=1,IDIM
  191. XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0*
  192. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  193. & XCOOR((NN3-1)*(IDIM+1)+I1))
  194. XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0*
  195. & (XCOOR((NN3-1)*(IDIM+1)+I1)+
  196. & XCOOR((NN5-1)*(IDIM+1)+I1))
  197. XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0*
  198. & (XCOOR((NN5-1)*(IDIM+1)+I1)+
  199. & XCOOR((NN1-1)*(IDIM+1)+I1))
  200. XCOOR((NN7-1)*(IDIM+1)+I1)=0.5D0*
  201. & (XCOOR((NN10-1)*(IDIM+1)+I1)+
  202. & XCOOR((NN1-1)*(IDIM+1)+I1))
  203. XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0*
  204. & (XCOOR((NN10-1)*(IDIM+1)+I1)+
  205. & XCOOR((NN3-1)*(IDIM+1)+I1))
  206. XCOOR((NN9-1)*(IDIM+1)+I1)=0.5D0*
  207. & (XCOOR((NN10-1)*(IDIM+1)+I1)+
  208. & XCOOR((NN5-1)*(IDIM+1)+I1))
  209. XCOOR((NN11-1)*(IDIM+1)+I1)=
  210. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  211. & XCOOR((NN3-1)*(IDIM+1)+I1)+
  212. & XCOOR((NN5-1)*(IDIM+1)+I1))/3.0D0
  213. XCOOR((NN12-1)*(IDIM+1)+I1)=
  214. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  215. & XCOOR((NN3-1)*(IDIM+1)+I1)+
  216. & XCOOR((NN10-1)*(IDIM+1)+I1))/3.0D0
  217. XCOOR((NN13-1)*(IDIM+1)+I1)=
  218. & (XCOOR((NN3-1)*(IDIM+1)+I1)+
  219. & XCOOR((NN5-1)*(IDIM+1)+I1)+
  220. & XCOOR((NN10-1)*(IDIM+1)+I1))/3.0D0
  221. XCOOR((NN14-1)*(IDIM+1)+I1)=
  222. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  223. & XCOOR((NN5-1)*(IDIM+1)+I1)+
  224. & XCOOR((NN10-1)*(IDIM+1)+I1))/3.0D0
  225. XCOOR((NN15-1)*(IDIM+1)+I1)=
  226. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  227. & XCOOR((NN3-1)*(IDIM+1)+I1)+
  228. & XCOOR((NN5-1)*(IDIM+1)+I1)+
  229. & XCOOR((NN10-1)*(IDIM+1)+I1))/4.0D0
  230. ENDDO
  231. ENDDO
  232. ELSEIF(IPT1.ITYPEL .EQ. 36)THEN
  233. C
  234. C********** PY19
  235. C
  236. IF(LOPY19)THEN
  237. C Elt already referred
  238. WRITE(IOIMP,*) 'Subroutine kdom2'
  239. WRITE(IOIMP,*) 'Mesh type not recognized'
  240. CALL ERREUR(5)
  241. ENDIF
  242. LOPY19=.TRUE.
  243. C
  244. DO IEL=1,NBEL,1
  245. NN1=IPT1.NUM(1,IEL)
  246. NN2=IPT1.NUM(2,IEL)
  247. NN3=IPT1.NUM(3,IEL)
  248. NN4=IPT1.NUM(4,IEL)
  249. NN5=IPT1.NUM(5,IEL)
  250. NN6=IPT1.NUM(6,IEL)
  251. NN7=IPT1.NUM(7,IEL)
  252. NN8=IPT1.NUM(8,IEL)
  253. NN9=IPT1.NUM(9,IEL)
  254. NN10=IPT1.NUM(10,IEL)
  255. NN11=IPT1.NUM(11,IEL)
  256. NN12=IPT1.NUM(12,IEL)
  257. NN13=IPT1.NUM(13,IEL)
  258. NN14=IPT1.NUM(14,IEL)
  259. NN15=IPT1.NUM(15,IEL)
  260. NN16=IPT1.NUM(16,IEL)
  261. NN17=IPT1.NUM(17,IEL)
  262. NN18=IPT1.NUM(18,IEL)
  263. NN19=IPT1.NUM(19,IEL)
  264. C Computation of the center position
  265. CALL KDOM7(NN1,NN3,NN5,NN7,XCEN)
  266. IF(IERR.NE.0)GOTO 9999
  267. DO I1=1,IDIM
  268. XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0*
  269. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  270. & XCOOR((NN3-1)*(IDIM+1)+I1))
  271. XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0*
  272. & (XCOOR((NN3-1)*(IDIM+1)+I1)+
  273. & XCOOR((NN5-1)*(IDIM+1)+I1))
  274. XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0*
  275. & (XCOOR((NN5-1)*(IDIM+1)+I1)+
  276. & XCOOR((NN7-1)*(IDIM+1)+I1))
  277. XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0*
  278. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  279. & XCOOR((NN7-1)*(IDIM+1)+I1))
  280. XCOOR((NN9-1)*(IDIM+1)+I1)=0.5D0*
  281. & (XCOOR((NN13-1)*(IDIM+1)+I1)+
  282. & XCOOR((NN1-1)*(IDIM+1)+I1))
  283. XCOOR((NN10-1)*(IDIM+1)+I1)=0.5D0*
  284. & (XCOOR((NN13-1)*(IDIM+1)+I1)+
  285. & XCOOR((NN3-1)*(IDIM+1)+I1))
  286. XCOOR((NN11-1)*(IDIM+1)+I1)=0.5D0*
  287. & (XCOOR((NN13-1)*(IDIM+1)+I1)+
  288. & XCOOR((NN5-1)*(IDIM+1)+I1))
  289. XCOOR((NN12-1)*(IDIM+1)+I1)=0.5D0*
  290. & (XCOOR((NN13-1)*(IDIM+1)+I1)+
  291. & XCOOR((NN7-1)*(IDIM+1)+I1))
  292. XCOOR((NN14-1)*(IDIM+1)+I1)=XCEN(I1)
  293. XCOOR((NN15-1)*(IDIM+1)+I1)=
  294. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  295. & XCOOR((NN3-1)*(IDIM+1)+I1)+
  296. & XCOOR((NN13-1)*(IDIM+1)+I1))/3.0D0
  297. XCOOR((NN16-1)*(IDIM+1)+I1)=
  298. & (XCOOR((NN3-1)*(IDIM+1)+I1)+
  299. & XCOOR((NN5-1)*(IDIM+1)+I1)+
  300. & XCOOR((NN13-1)*(IDIM+1)+I1))/3.0D0
  301. XCOOR((NN17-1)*(IDIM+1)+I1)=
  302. & (XCOOR((NN5-1)*(IDIM+1)+I1)+
  303. & XCOOR((NN7-1)*(IDIM+1)+I1)+
  304. & XCOOR((NN13-1)*(IDIM+1)+I1))/3.0D0
  305. XCOOR((NN18-1)*(IDIM+1)+I1)=
  306. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  307. & XCOOR((NN7-1)*(IDIM+1)+I1)+
  308. & XCOOR((NN13-1)*(IDIM+1)+I1))/3.0D0
  309. ENDDO
  310. C
  311. C************* To compute the center, we divide the pyramid into
  312. C 4 tetrahedra (see kdom3)
  313. C
  314. CALL KDOM3(NN1,NN3,NN5,NN7,NN14,NN13,XCEN,VOL)
  315. DO I1=1,3,1
  316. XCOOR((NN19-1)*(IDIM+1)+I1)=XCEN(I1)
  317. ENDDO
  318. ENDDO
  319. ELSEIF(IPT1.ITYPEL .EQ. 34)THEN
  320. C
  321. C********** PR21
  322. C
  323. IF(LOPR21)THEN
  324. C Elt already referred
  325. WRITE(IOIMP,*) 'Subroutine kdom2'
  326. WRITE(IOIMP,*) 'Mesh type not recognized'
  327. CALL ERREUR(5)
  328. ENDIF
  329. LOPR21=.TRUE.
  330. C
  331. DO IEL=1,NBEL,1
  332. NN1=IPT1.NUM(1,IEL)
  333. NN2=IPT1.NUM(2,IEL)
  334. NN3=IPT1.NUM(3,IEL)
  335. NN4=IPT1.NUM(4,IEL)
  336. NN5=IPT1.NUM(5,IEL)
  337. NN6=IPT1.NUM(6,IEL)
  338. NN7=IPT1.NUM(7,IEL)
  339. NN8=IPT1.NUM(8,IEL)
  340. NN9=IPT1.NUM(9,IEL)
  341. NN10=IPT1.NUM(10,IEL)
  342. NN11=IPT1.NUM(11,IEL)
  343. NN12=IPT1.NUM(12,IEL)
  344. NN13=IPT1.NUM(13,IEL)
  345. NN14=IPT1.NUM(14,IEL)
  346. NN15=IPT1.NUM(15,IEL)
  347. NN16=IPT1.NUM(16,IEL)
  348. NN17=IPT1.NUM(17,IEL)
  349. NN18=IPT1.NUM(18,IEL)
  350. NN19=IPT1.NUM(19,IEL)
  351. NN20=IPT1.NUM(20,IEL)
  352. NN21=IPT1.NUM(21,IEL)
  353. DO I1=1,IDIM
  354. XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0*
  355. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  356. & XCOOR((NN3-1)*(IDIM+1)+I1))
  357. XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0*
  358. & (XCOOR((NN3-1)*(IDIM+1)+I1)+
  359. & XCOOR((NN5-1)*(IDIM+1)+I1))
  360. XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0*
  361. & (XCOOR((NN5-1)*(IDIM+1)+I1)+
  362. & XCOOR((NN1-1)*(IDIM+1)+I1))
  363. XCOOR((NN7-1)*(IDIM+1)+I1)=0.5D0*
  364. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  365. & XCOOR((NN10-1)*(IDIM+1)+I1))
  366. XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0*
  367. & (XCOOR((NN3-1)*(IDIM+1)+I1)+
  368. & XCOOR((NN12-1)*(IDIM+1)+I1))
  369. XCOOR((NN9-1)*(IDIM+1)+I1)=0.5D0*
  370. & (XCOOR((NN5-1)*(IDIM+1)+I1)+
  371. & XCOOR((NN14-1)*(IDIM+1)+I1))
  372. XCOOR((NN11-1)*(IDIM+1)+I1)=0.5D0*
  373. & (XCOOR((NN10-1)*(IDIM+1)+I1)+
  374. & XCOOR((NN12-1)*(IDIM+1)+I1))
  375. XCOOR((NN13-1)*(IDIM+1)+I1)=0.5D0*
  376. & (XCOOR((NN12-1)*(IDIM+1)+I1)+
  377. & XCOOR((NN14-1)*(IDIM+1)+I1))
  378. XCOOR((NN15-1)*(IDIM+1)+I1)=0.5D0*
  379. & (XCOOR((NN10-1)*(IDIM+1)+I1)+
  380. & XCOOR((NN14-1)*(IDIM+1)+I1))
  381. XCOOR((NN19-1)*(IDIM+1)+I1)=
  382. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  383. & XCOOR((NN3-1)*(IDIM+1)+I1)+
  384. & XCOOR((NN5-1)*(IDIM+1)+I1))/3.0D0
  385. XCOOR((NN20-1)*(IDIM+1)+I1)=
  386. & (XCOOR((NN10-1)*(IDIM+1)+I1)+
  387. & XCOOR((NN12-1)*(IDIM+1)+I1)+
  388. & XCOOR((NN14-1)*(IDIM+1)+I1))/3.0D0
  389. ENDDO
  390. C Computation of the QUA center positions
  391. CALL KDOM7(NN1,NN10,NN12,NN3,XCEN)
  392. IF(IERR.NE.0)GOTO 9999
  393. DO I1=1,IDIM,1
  394. XCOOR((NN16-1)*(IDIM+1)+I1)=XCEN(I1)
  395. ENDDO
  396. CALL KDOM7(NN5,NN3,NN12,NN14,XCEN)
  397. IF(IERR.NE.0)GOTO 9999
  398. DO I1=1,IDIM,1
  399. XCOOR((NN17-1)*(IDIM+1)+I1)=XCEN(I1)
  400. ENDDO
  401. CALL KDOM7(NN10,NN1,NN5,NN14,XCEN)
  402. IF(IERR.NE.0)GOTO 9999
  403. DO I1=1,IDIM,1
  404. XCOOR((NN18-1)*(IDIM+1)+I1)=XCEN(I1)
  405. ENDDO
  406. C
  407. C************* To compute the center, we divide the PR21
  408. C into 2 tetrahedra and 3 pyramids
  409. C (1,3,5,PC), (10,14,12,PC), (1,10,12,3,16,PC),
  410. C (3,12,14,5,17,PC),(1,5,14,10,18,PC)
  411. C |
  412. C center of 3,12,14,5
  413. C The results does not depend of the choice of
  414. C PC.
  415. C
  416. C
  417. C
  418. DO I1=1,3,1
  419. XCOOR((NN21-1)*(IDIM+1)+I1)=
  420. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  421. & XCOOR((NN3-1)*(IDIM+1)+I1)+
  422. & XCOOR((NN5-1)*(IDIM+1)+I1)+
  423. & XCOOR((NN10-1)*(IDIM+1)+I1)+
  424. & XCOOR((NN12-1)*(IDIM+1)+I1)+
  425. & XCOOR((NN14-1)*(IDIM+1)+I1))/6.0D0
  426.  
  427. ENDDO
  428. C
  429. CALL KDOM3(NN10,NN12,NN3,NN1,NN16,NN21,X1,VOL1)
  430. CALL KDOM3(NN3,NN12,NN14,NN5,NN17,NN21,X2,VOL2)
  431. CALL KDOM3(NN1,NN5,NN14,NN10,NN18,NN21,X3,VOL3)
  432. CALL KDOM4(NN1,NN3,NN5,NN21,X4,VOL4)
  433. CALL KDOM4(NN10,NN14,NN12,NN21,X5,VOL5)
  434. VOL=VOL1+VOL2+VOL3+VOL4+VOL5
  435. C
  436. DO I1=1,3,1
  437. XCOOR((NN21-1)*(IDIM+1)+I1)=(VOL1*X1(I1)+VOL2*X2(I1)
  438. & +VOL3*X3(I1)+VOL4*X4(I1)+VOL5*X5(I1))/VOL
  439. ENDDO
  440. C
  441. C********** We check whether the center of gravity of the PR21 is
  442. C inside its domain
  443. C
  444. CALL KDOM8(NN1,NN2,NN3,NN4,NN5,NN6,NN7,NN8,NN9,NN10,
  445. & NN11,NN12,NN13,NN14,NN15,NN16,NN17,NN18,NN19,NN20,
  446. & NN21)
  447. IF(IERR.NE.0)GOTO 9999
  448. ENDDO
  449. IF(NBSOUS .NE. 1) SEGDES IPT1
  450. ELSEIF(IPT1.ITYPEL .EQ. 33)THEN
  451. C
  452. C********** CU27
  453. C
  454. IF(LOCU27)THEN
  455. C Elt already referred
  456. WRITE(IOIMP,*) 'Subroutine kdom2'
  457. WRITE(IOIMP,*) 'Mesh type not recognized'
  458. CALL ERREUR(5)
  459. ENDIF
  460. LOCU27=.TRUE.
  461. C
  462. DO IEL=1,NBEL,1
  463. NN1=IPT1.NUM(1,IEL)
  464. NN2=IPT1.NUM(2,IEL)
  465. NN3=IPT1.NUM(3,IEL)
  466. NN4=IPT1.NUM(4,IEL)
  467. NN5=IPT1.NUM(5,IEL)
  468. NN6=IPT1.NUM(6,IEL)
  469. NN7=IPT1.NUM(7,IEL)
  470. NN8=IPT1.NUM(8,IEL)
  471. NN9=IPT1.NUM(9,IEL)
  472. NN10=IPT1.NUM(10,IEL)
  473. NN11=IPT1.NUM(11,IEL)
  474. NN12=IPT1.NUM(12,IEL)
  475. NN13=IPT1.NUM(13,IEL)
  476. NN14=IPT1.NUM(14,IEL)
  477. NN15=IPT1.NUM(15,IEL)
  478. NN16=IPT1.NUM(16,IEL)
  479. NN17=IPT1.NUM(17,IEL)
  480. NN18=IPT1.NUM(18,IEL)
  481. NN19=IPT1.NUM(19,IEL)
  482. NN20=IPT1.NUM(20,IEL)
  483. NN21=IPT1.NUM(21,IEL)
  484. NN22=IPT1.NUM(22,IEL)
  485. NN23=IPT1.NUM(23,IEL)
  486. NN24=IPT1.NUM(24,IEL)
  487. NN25=IPT1.NUM(25,IEL)
  488. NN26=IPT1.NUM(26,IEL)
  489. NN27=IPT1.NUM(27,IEL)
  490. DO I1=1,IDIM
  491. XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0*
  492. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  493. & XCOOR((NN3-1)*(IDIM+1)+I1))
  494. XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0*
  495. & (XCOOR((NN3-1)*(IDIM+1)+I1)+
  496. & XCOOR((NN5-1)*(IDIM+1)+I1))
  497. XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0*
  498. & (XCOOR((NN5-1)*(IDIM+1)+I1)+
  499. & XCOOR((NN7-1)*(IDIM+1)+I1))
  500. XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0*
  501. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  502. & XCOOR((NN7-1)*(IDIM+1)+I1))
  503. XCOOR((NN9-1)*(IDIM+1)+I1)=0.5D0*
  504. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  505. & XCOOR((NN13-1)*(IDIM+1)+I1))
  506. XCOOR((NN10-1)*(IDIM+1)+I1)=0.5D0*
  507. & (XCOOR((NN3-1)*(IDIM+1)+I1)+
  508. & XCOOR((NN15-1)*(IDIM+1)+I1))
  509. XCOOR((NN11-1)*(IDIM+1)+I1)=0.5D0*
  510. & (XCOOR((NN5-1)*(IDIM+1)+I1)+
  511. & XCOOR((NN17-1)*(IDIM+1)+I1))
  512. XCOOR((NN12-1)*(IDIM+1)+I1)=0.5D0*
  513. & (XCOOR((NN7-1)*(IDIM+1)+I1)+
  514. & XCOOR((NN19-1)*(IDIM+1)+I1))
  515. XCOOR((NN14-1)*(IDIM+1)+I1)=0.5D0*
  516. & (XCOOR((NN13-1)*(IDIM+1)+I1)+
  517. & XCOOR((NN15-1)*(IDIM+1)+I1))
  518. XCOOR((NN16-1)*(IDIM+1)+I1)=0.5D0*
  519. & (XCOOR((NN17-1)*(IDIM+1)+I1)+
  520. & XCOOR((NN15-1)*(IDIM+1)+I1))
  521. XCOOR((NN18-1)*(IDIM+1)+I1)=0.5D0*
  522. & (XCOOR((NN17-1)*(IDIM+1)+I1)+
  523. & XCOOR((NN19-1)*(IDIM+1)+I1))
  524. XCOOR((NN20-1)*(IDIM+1)+I1)=0.5D0*
  525. & (XCOOR((NN19-1)*(IDIM+1)+I1)+
  526. & XCOOR((NN13-1)*(IDIM+1)+I1))
  527. ENDDO
  528. C Computation of the QUA center positions
  529. CALL KDOM7(NN1,NN13,NN15,NN3,XCEN)
  530. IF(IERR.NE.0) GOTO 9999
  531. DO I1=1,IDIM,1
  532. XCOOR((NN21-1)*(IDIM+1)+I1)=XCEN(I1)
  533. ENDDO
  534. CALL KDOM7(NN3,NN15,NN17,NN5,XCEN)
  535. IF(IERR.NE.0)GOTO 9999
  536. DO I1=1,IDIM,1
  537. XCOOR((NN22-1)*(IDIM+1)+I1)=XCEN(I1)
  538. ENDDO
  539. CALL KDOM7(NN5,NN17,NN19,NN7,XCEN)
  540. IF(IERR.NE.0)GOTO 9999
  541. DO I1=1,IDIM,1
  542. XCOOR((NN23-1)*(IDIM+1)+I1)=XCEN(I1)
  543. ENDDO
  544. CALL KDOM7(NN7,NN19,NN13,NN1,XCEN)
  545. IF(IERR.NE.0)GOTO 9999
  546. DO I1=1,IDIM,1
  547. XCOOR((NN24-1)*(IDIM+1)+I1)=XCEN(I1)
  548. ENDDO
  549. CALL KDOM7(NN1,NN3,NN5,NN7,XCEN)
  550. IF(IERR.NE.0)GOTO 9999
  551. DO I1=1,IDIM,1
  552. XCOOR((NN25-1)*(IDIM+1)+I1)=XCEN(I1)
  553. ENDDO
  554. CALL KDOM7(NN13,NN19,NN17,NN15,XCEN)
  555. IF(IERR.NE.0)GOTO 9999
  556. DO I1=1,IDIM,1
  557. XCOOR((NN26-1)*(IDIM+1)+I1)=XCEN(I1)
  558. ENDDO
  559. C
  560. C************* To compute the center, we divide the cube
  561. C into 6 pyramids.
  562. C (1,3,5,7,25,PC)
  563. C (13,19,17,15,26,PC)
  564. C (1,13,15,3,21,PC)
  565. C (3,15,17,5,22,PC)
  566. C (5,17,19,7,23,PC)
  567. C (1,7,19,13,24,PC)
  568.  
  569. C The results does not depend of the choice of
  570. C PC.
  571. C
  572. DO I1=1,3,1
  573. XCOOR((NN27-1)*(IDIM+1)+I1)=
  574. & (XCOOR((NN1-1)*(IDIM+1)+I1)+
  575. & XCOOR((NN3-1)*(IDIM+1)+I1)+
  576. & XCOOR((NN5-1)*(IDIM+1)+I1)+
  577. & XCOOR((NN7-1)*(IDIM+1)+I1)+
  578. & XCOOR((NN13-1)*(IDIM+1)+I1)+
  579. & XCOOR((NN15-1)*(IDIM+1)+I1)+
  580. & XCOOR((NN17-1)*(IDIM+1)+I1)+
  581. & XCOOR((NN19-1)*(IDIM+1)+I1))/8.0D0
  582.  
  583. ENDDO
  584. C
  585. CALL KDOM3(NN1,NN3,NN5,NN7,NN25,NN27,X1,VOL1)
  586. CALL KDOM3(NN13,NN19,NN17,NN15,NN26,NN27,X2,VOL2)
  587. CALL KDOM3(NN1,NN13,NN15,NN3,NN21,NN27,X3,VOL3)
  588. CALL KDOM3(NN3,NN15,NN17,NN5,NN22,NN27,X4,VOL4)
  589. CALL KDOM3(NN5,NN17,NN19,NN7,NN23,NN27,X5,VOL5)
  590. CALL KDOM3(NN1,NN7,NN19,NN13,NN24,NN27,X6,VOL6)
  591. VOL=VOL1+VOL2+VOL3+VOL4+VOL5+VOL6
  592. DO I1=1,3,1
  593. XCOOR((NN27-1)*(IDIM+1)+I1)=(VOL1*X1(I1)+VOL2*X2(I1)
  594. & +VOL3*X3(I1)+VOL4*X4(I1)+VOL5*X5(I1)+VOL6*X6(I1))
  595. & /VOL
  596. ENDDO
  597. C
  598. C********** We check whether the center of gravity of the CU27 is
  599. C inside its domain
  600. C
  601. CALL KDOM9(NN1,NN2,NN3,NN4,NN5,NN6,NN7,NN8,NN9,NN10,
  602. & NN11,NN12,NN13,NN14,NN15,NN16,NN17,NN18,NN19,NN20,
  603. & NN21,NN22,NN23,NN24,NN25,NN26,NN27)
  604. IF(IERR.NE.0)GOTO 9999
  605.  
  606. ENDDO
  607. ELSE
  608. WRITE(IOIMP,*) 'Euler model'
  609. CALL ERREUR(21)
  610. ENDIF
  611. IF(NBSOUS .NE. 1) SEGDES IPT1
  612. ENDDO
  613. C
  614. SEGDES MELEME
  615. C
  616. 9999 RETURN
  617. C
  618. END
  619.  
  620.  
  621.  
  622.  

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