Télécharger hhomat.eso

Retour à la liste

Numérotation des lignes :

hhomat
  1. C HHOMAT SOURCE OF166741 24/06/19 21:15:07 11942
  2. C HHOMAT SOURCE
  3.  
  4. SUBROUTINE HHOMAT (IPMODL, IPCARA, iret)
  5.  
  6. IMPLICIT INTEGER(I-N)
  7. IMPLICIT REAL*8(A-H,O-Z)
  8.  
  9. -INC PPARAM
  10. -INC CCOPTIO
  11.  
  12. -INC CCHHOPA
  13. -INC CCHHOPR
  14.  
  15. -INC SMCOORD
  16. -INC SMELEME
  17. -INC SMMODEL
  18. -INC SMCHAML
  19. -INC SMLENTI
  20. -INC SMLREEL
  21.  
  22. EXTERNAL LONG
  23.  
  24. CHARACTER*(LOCHAI) chaHHO, charerr
  25. CHARACTER*(NCONCH) conHHO
  26. CHARACTER*(LOCOMP) mocomp
  27.  
  28. iret = 0
  29.  
  30. mmodel = IPMODL
  31. c* segact,mmodel*nomod (segment actif en entree)
  32. NSOUM = mmodel.kmodel(/1)
  33.  
  34. mchelm = IPCARA
  35. N1 = mchelm.IMACHE(/1)
  36. c* segact,mchelm
  37.  
  38. DO im = 1, NSOUM
  39. imodel = mmodel.kmodel(im)
  40.  
  41. IF (imodel.nefmod .NE. HHO_NUM_ELEMENT) GOTO 100
  42. CALL HHONOB(imodel, nobHHO, iret)
  43. IF (nobHHO.LE.0) THEN
  44. write(ioimp,*) 'HHOMAT: nobHHO undefined'
  45. CALL ERREUR(5)
  46. RETURN
  47. END IF
  48.  
  49. iva = imodel.IVAMOD(nobHHO)
  50. CALL QUEVAL(iva,'MOT ',iret,lchho,r_z,chaHHO,b_z,i_z)
  51. if (iret .NE. 0) CALL ERREUR(5)
  52.  
  53. imaHHO = imodel.IMAMOD
  54. CALL PLACE2(mchelm.IMACHE(1),N1,i_z,imaHHO)
  55. IF (i_z .EQ. 0) GOTO 100
  56. conHHO = imodel.CONMOD
  57. IF (conHHO.NE.mchelm.CONCHE(i_z)) GOTO 100
  58.  
  59. ISUPM = mchelm.INFCHE(i_z,6)
  60. IF (ISUPM.NE.3) THEN
  61. write(ioimp,*) 'HHOMAT: ISUPM /= 3 (',ISUPM,')'
  62. END IF
  63. IF (mchelm.INFCHE(i_z,4).NE.imodel.infmod(2+ISUPM)) THEN
  64. write(ioimp,*) 'HHOMAT: Pointeur infche/infmod ???'
  65. END IF
  66.  
  67. mchaml = mchelm.ICHAML(i_z)
  68. c* segact,mchaml
  69. noCOMP = mchaml.IELVAL(/1)
  70. moCOMP = 'STAB '
  71. CALL PLACE(mchaml.NOMCHE(1),noCOMP,i_z,moCOMP)
  72. c* pas de composante STAB, je plante tout !
  73. IF (i_z .EQ. 0) THEN
  74. write(ioimp,*) 'STAB component is needed for HHO formulation'
  75. CALL ERREUR(21)
  76. GOTO 100
  77. END IF
  78.  
  79. meleme = imaHHO
  80. c* segact,meleme
  81. NBNOE = meleme.num(/1)
  82. NBELT = meleme.num(/2)
  83.  
  84. mlenti = imodel.ivamod(nobhho+1)
  85. c* segact,mlenti
  86.  
  87. NBFAC = mlenti.lect(7)
  88. c- Le nombre de points d'integration correspond a ISUPM !
  89. NBPIN = imodel.infele(6)
  90.  
  91. mlent3 = imodel.ivamod(nobhho+3)
  92. c* segact,mlent3
  93. nb3 = mlent3.lect(/1) / 2
  94. nbfac3 = 2 * NBFAC
  95.  
  96. c-dbg write(ioimp,*) 'HHOMAT',im,NSOUM,nbnoe,nbelt,nbpin
  97. IF (IDIM .NE. mlenti.lect(1)) then
  98. write(ioimp,*) 'HHOMAT: Bizarre idim'
  99. END IF
  100. IF (NBNOE .NE. mlenti.lect(6)) then
  101. write(ioimp,*) 'HHOMAT: Bizarre nbnoe'
  102. END IF
  103. IF (IDIM.EQ.2 .AND. NBFAC .NE. NBNOE) THEN
  104. write(ioimp,*) 'HHOMAT: Bizarre nbfac2d'
  105. END IF
  106. IF (NBPIN .NE. mlenti.lect(8)) then
  107. write(ioimp,*) 'HHOMAT: Bizarre nb.pint'
  108. END IF
  109. IF (nb3.ne.(NBELT*NBFAC)) THEN
  110. write(ioimp,*) 'HHOMAT: Bizarre mlent3'
  111. END IF
  112.  
  113. N2 = noCOMP
  114. moCOMP = 'PIHO '
  115. CALL PLACE(mchaml.NOMCHE(1),noCOMP,icpi,moCOMP)
  116. IF (icpi .EQ. 0) THEN
  117. N2 = N2 + 1
  118. icpi = N2
  119. END IF
  120. moCOMP = 'BHHO '
  121. CALL PLACE(mchaml.NOMCHE(1),noCOMP,icmb,moCOMP)
  122. IF (icmb .EQ. 0) THEN
  123. N2 = N2 + 1
  124. icmb = N2
  125. END IF
  126. moCOMP = 'SHHO '
  127. CALL PLACE(mchaml.NOMCHE(1),noCOMP,icms,moCOMP)
  128. IF (icms .EQ. 0) THEN
  129. N2 = N2 + 1
  130. icms = N2
  131. END IF
  132. IF (N2.GT.noCOMP) THEN
  133. segadj,mchaml
  134. ELSE
  135. segact,mchaml*MOD
  136. END IF
  137.  
  138. mchaml.NOMCHE(icpi) = 'PIHO '
  139. mchaml.TYPCHE(icpi) = 'REAL*8 '
  140. N1PTEL = NBPIN
  141. N1EL = NBELT
  142. N2PTEL = 0
  143. N2EL = 0
  144. IF (icpi.GT.nocomp) THEN
  145. SEGINI,melval
  146. mchaml.IELVAL(icpi) = melval
  147. ELSE
  148. melval = mchaml.IELVAL(icpi)
  149. SEGADJ,melval
  150. END IF
  151. IVPIHO = melval
  152.  
  153. mchaml.NOMCHE(icmb) = 'BHHO '
  154. mchaml.TYPCHE(icmb) = 'POINTEURLISTREEL'
  155. N1PTEL = 0
  156. N1EL = 0
  157. N2PTEL = NBPIN
  158. N2EL = NBELT
  159. IF (icmb.GT.nocomp) THEN
  160. SEGINI,melval
  161. mchaml.IELVAL(icmb) = melval
  162. ELSE
  163. melval = mchaml.IELVAL(icmb)
  164. SEGADJ,melval
  165. END IF
  166. IVMBHO = melval
  167.  
  168. mchaml.NOMCHE(icms) = 'SHHO '
  169. mchaml.TYPCHE(icms) = 'POINTEURLISTREEL'
  170. melval = mchaml.IELVAL(icms)
  171. N1PTEL = 0
  172. N1EL = 0
  173. N2PTEL = 1
  174. N2EL = NBELT
  175. IF (icms.GT.nocomp) THEN
  176. SEGINI,melval
  177. mchaml.IELVAL(icms) = melval
  178. ELSE
  179. melval = mchaml.IELVAL(icms)
  180. SEGADJ,melval
  181. END IF
  182. IVMSTA = melval
  183.  
  184. C- Preparation des tableaux pour l'appel a HHOC3M et les fonctions de
  185. C- calcul associees a l'element HHO ici present
  186. ne_0 = 9
  187. * connectivity des faces : en 2D 2 noeuds par aretes --> en 3D a faire !
  188. ne_co = nbfac3
  189. ie_0 = 0
  190. ie_co = ie_0 + ne_0
  191. ie_1 = ie_co + ne_co
  192.  
  193. nr_0 = 0
  194. nr_co = IDIM * nbnoe
  195. nr_pi = NBPIN
  196. nr_bh = mlenti.lect(14)
  197. nr_sh = mlenti.lect(16)
  198. ir_0 = 0
  199. ir_co = ir_0 + nr_0
  200. ir_1 = ir_co + nr_co
  201.  
  202. JG = ie_1
  203. SEGINI,mlent2
  204. ile = JG
  205.  
  206. mlent2.lect(ie_0 + 1) = ie_co + 1
  207. mlent2.lect(ie_0 + 2) = ie_1
  208. mlent2.lect(ie_0 + 3) = ir_co + 1
  209. mlent2.lect(ie_0 + 4) = ir_1
  210.  
  211. JG = nr_0 + nr_co + MAX(nr_pi,nr_bh,nr_sh)
  212. SEGINI,mlree2
  213. ilr = JG
  214.  
  215. c* SEGACT,MCOORD*NOMOD
  216.  
  217. DO IEL = 1, NBELT
  218.  
  219. C- Recuperation des coordonnees des noeuds de l element IEL
  220. CALL HHOCOO(meleme.num,NBNOE,mcoord.xcoor(1), IEL,
  221. & mlree2.prog(ir_co+1),iret)
  222. IF (iret.NE.0) RETURN
  223. c-dbg write(ioimp,*) 'HHOMAT-HHOCOO',IEL,NBNOE,ir_co,ir_1
  224. c-dbg write(ioimp,*) ' ',(meleme.num(i,IEL),i=1,NBNOE)
  225. c-dbg write(ioimp,*) ' X',(mlree2.prog(ir_co+i),i=1,NBNOE)
  226. c-dbg write(ioimp,*) ' Y',(mlree2.prog(ir_co+NBNOE+i),i=1,NBNOE)
  227. C- Connectivite des faces de l element IELT en tenant compte du signe :
  228. C- il faudra faire une fonction pour le 3D. pour le 2D que des SEG2
  229. C- Dans la bibliotheque, la numerotation locale des sommets d'une cellule
  230. C- debute a 0 ! Attention en 2D sur le dernier element !
  231. in1 = (IEL-1) * NBFAC3
  232. DO i = 1, NBFAC
  233. j = 2 * i
  234. in2 = in1 + j
  235. if (mlent3.lect(in2-1).ne.2)
  236. & write(ioimp,*) 'HHOMAT : Bizarre mlent3 2',iel,i,in1,in2
  237. ie2 = mlent3.lect(in2)
  238. if (ie2.eq.0)
  239. & write(ioimp,*) 'HHOMAT : Bizarre mlent3 ie2',iel,i,in1,in2
  240. IF (ie2.GT.0) THEN
  241. mlent2.lect(ie_co+j-1) = i-1
  242. mlent2.lect(ie_co+j ) = MOD(i,NBFAC)
  243. ELSE
  244. mlent2.lect(ie_co+j-1) = MOD(i,NBFAC)
  245. mlent2.lect(ie_co+j ) = i-1
  246. END IF
  247. c-dbg write(ioimp,*) 'hhomat - iel',iel,i,in2,ie2
  248. END DO
  249.  
  250. C- Recuperation des poids de Gauss :
  251. mlent2.lect(ie_0 + 5) = nr_pi
  252. mlent2.lect(ie_0 + 6) = ir_1
  253. mlent2.lect(ie_0 + 7) = 0
  254.  
  255. DO i = ir_1+1, ilr
  256. mlree2.prog(i) = 0.D0
  257. END DO
  258.  
  259. iretc = 0
  260. CALL HHOC3M('INTG',chaHHO(1:lchho),
  261. & HHO_NomLib, HHO_MaxLib,
  262. & mlent2.lect,ile, mlree2.prog,ilr,
  263. & iretc,charerr)
  264. IF (iretc.NE.0) THEN
  265. write(ioimp,*) 'HHOMAT -> HHOC3M -> INTG : ERROR ='
  266. write(ioimp,*) charerr(1:LONG(charerr))
  267. iret = 21
  268. return
  269. END IF
  270.  
  271. melval = IVPIHO
  272. DO IGA = 1, NBPIN
  273. melval.velche(IGA,IEL) = mlree2.prog(ir_1+iga)
  274. END DO
  275.  
  276. C- Recuperation des matrices B en chaque point de Gauss :
  277. mlent2.lect(ie_0 + 5) = ir_1 + 1
  278. mlent2.lect(ie_0 + 6) = ir_1 + nr_bh
  279.  
  280. melval = IVMBHO
  281. DO IGA = 1, NBPIN
  282.  
  283. mlent2.lect(ie_0 + 7) = iga
  284. DO i = ir_1+1, ilr
  285. mlree2.prog(i) = 0.D0
  286. END DO
  287.  
  288. iretc = 0
  289. CALL HHOC3M('BHHO',chaHHO(1:lchho),
  290. & HHO_NomLib, HHO_MaxLib,
  291. & mlent2.lect,ile, mlree2.prog,ilr,
  292. & iretc,charerr)
  293. IF (iretc.NE.0) THEN
  294. write(ioimp,*) 'HHOMAT -> HHOC3M -> BHHO : ERROR ='
  295. write(ioimp,*) charerr(1:LONG(charerr))
  296. iret = 21
  297. return
  298. END IF
  299.  
  300. c- La matrice BHHO(9,NBDDL) :
  301. NLIG = 9
  302. NCOL = mlenti.lect(11)
  303. JG = nr_bh
  304. SEGINI,mlreel
  305. if (JG.ne.(NLIG*NCOL)) write(ioimp,*)'HHOMAT - Bizarre nr_bh',iga
  306. c* Matrice BHHO stockee ligne par ligne dans la bibliotheque.
  307. c* Matrice BHHO stockee colonne par colonne dans Cast3M.
  308. DO j = 1, NCOL
  309. jc = NLIG * (j-1)
  310. DO i = 1, NLIG
  311. ic = NCOL * (i-1)
  312. mlreel.prog(jc+i) = mlree2.prog(ir_1+ic+j)
  313. END DO
  314. END DO
  315. c* SEGDES,mlreel
  316.  
  317. melval.ielche(IGA,IEL) = mlreel
  318.  
  319. END DO
  320.  
  321. C- Recuperation de la matrice de stabilisation :
  322. mlent2.lect(ie_0 + 5) = ir_1 + 1
  323. mlent2.lect(ie_0 + 6) = ir_1 + nr_sh
  324.  
  325. DO i = ir_1+1, ilr
  326. mlree2.prog(i) = 0.D0
  327. END DO
  328.  
  329. iretc = 0
  330. CALL HHOC3M('SHHO',chaHHO(1:lchho),
  331. & HHO_NomLib, HHO_MaxLib,
  332. & mlent2.lect,ile, mlree2.prog,ilr,
  333. & iretc,charerr)
  334. IF (iretc.NE.0) THEN
  335. write(ioimp,*) 'HHOMAT -> HHOC3M -> SHHO : ERROR ='
  336. write(ioimp,*) charerr(1:LONG(charerr))
  337. iret = 21
  338. return
  339. END IF
  340.  
  341. NLIG = mlenti.lect(11)
  342. NCOL = mlenti.lect(11)
  343. JG = nr_sh
  344. SEGINI,mlreel
  345. if (JG.ne.(NLIG*NCOL)) write(ioimp,*)'HHOMAT - Bizarre nr_sh'
  346. c* !! matrice SHHO stockee ligne par ligne dans la bibliotheque.
  347. c* !! matrice SHHO stockee colonne par colonne dans Cast3M.
  348. DO j = 1, NCOL
  349. jc = NLIG * (j-1)
  350. DO i = 1, NLIG
  351. ic = NCOL * (i-1)
  352. mlreel.prog(jc+i) = mlree2.prog(ir_1+ic+j)
  353. END DO
  354. END DO
  355. c* SEGDES,mlreel
  356.  
  357. melval = IVMSTA
  358. melval.ielche(1,IEL) = mlreel
  359.  
  360. END DO
  361.  
  362. 100 CONTINUE
  363. END DO
  364.  
  365. c RETURN
  366. END
  367.  
  368.  
  369.  
  370.  

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