Télécharger knrf.eso

Retour à la liste

Numérotation des lignes :

  1. C KNRF SOURCE BP208322 16/11/18 21:18:18 9177
  2. SUBROUTINE KNRF(MTABD,MCHELM,MCHPOI,MCHPO2)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. C
  6. C*************************************************************************
  7. C
  8. C OBJET : Cree deux CHAMPOINT faces et un mchaml
  9. C CHPC1 contenant la longueur des faces
  10. C ( faux dans le cas 3d )
  11. C CHPC2 contenant la normale aux faces.
  12. C MCHAML orientation (normale entrante ou sortante)
  13. C SYNTAXE : CHPC1 CHPC2 MCHAML = KNRF OBJDOM ;
  14. C
  15. C OBJDOM : TABLE de SOUSTYPE DOMAINE
  16. C
  17. C REMARQUE :
  18. C On verifie que la normale est dans le meme sens que pt1->pt3
  19. C des elements FACEL
  20. C Si pt1=pt3 on oriente la normale vers l'exterieur cad pt1->pt2
  21. C
  22. C Non teste si axisymetrie
  23. C*************************************************************************
  24. C MODIFICATION du 20/11/98
  25. C
  26. C Le calcul de la normale s'effectue a l'aide du jacobien calcule
  27. C par CALJBR (au lieu de CALJQB)
  28. C
  29. C*************************************************************************
  30. -INC CCOPTIO
  31. -INC CCGEOME
  32. -INC SMTABLE
  33. POINTEUR MTABD.MTABLE
  34. -INC SMELEME
  35. POINTEUR MELEMC.MELEME,MELEMD.MELEME,MELEMT.MELEME
  36. POINTEUR MELEMF.MELEME
  37. -INC SMCOORD
  38. -INC SMCHPOI
  39. -INC SMCHAML
  40. -INC SMLENTI
  41. -INC SIZFFB
  42. -INC CCREEL
  43. C
  44. PARAMETER (NBO=5)
  45. REAL*8 A2J(2,2,1),A3J(3,3,1),CFT,NORMX,NORMY,NORMZ
  46. REAL*8 LGR,LGR1,LGR2
  47. CHARACTER*8 TYPE,TYPC,NOM0,LIST1(NBO),LIST2(NBO)
  48. C Les elements complets ne sont pas tous implementes
  49. C On se ramene aux elements facep sans le centre de la face
  50. DATA LIST1/'SEG3 ','TRI4 ','TRI7 ','QUA5 ',
  51. & 'QUA9 '/
  52. DATA LIST2/'SEG2 ','TRI3 ','TRI6 ','QUA4 ',
  53. & 'QUA8 '/
  54. C
  55. C Acquisition des connectivites FACEL
  56. C write(6,*)' KNRF '
  57. TYPE=' '
  58. CALL ACMO(MTABD,'FACEL',TYPE,MELEMD)
  59. C write(6,*)' FACEL type=',type
  60. IF(TYPE.NE.'MAILLAGE')GO TO 90
  61. SEGACT MELEMD
  62. IPT2=MELEMD
  63. C
  64. C Acquisition des connectivites FACE
  65. C
  66. TYPE=' '
  67. CALL ACMO(MTABD,'FACE',TYPE,MELEMF)
  68. IF(TYPE.NE.'MAILLAGE')GO TO 90
  69. CALL KRIPAD(MELEMF,MLENTI)
  70. SEGACT MELEMF
  71. C
  72. C Creation du CHPOIN a 3 comp. pour la normale aux faces
  73. C
  74. NC=IDIM
  75. TYPE='FACE'
  76. CALL CRCHPT(TYPE,MELEMF,NC,MCHPOI)
  77. CALL LICHT(MCHPOI,MPOVAL,TYPC,IGEOM)
  78. C
  79. C Creation du CHPOIN pour la longueur
  80. C
  81. NC=1
  82. TYPE='FACE'
  83. CALL CRCHPT(TYPE,MELEMF,NC,MCHPO2)
  84. CALL LICHT(MCHPO2,MPOVA2,TYPC,IGEOM)
  85. SEGDES MELEMF
  86. C
  87. IAXI = 0
  88. IF(IFOMOD.EQ.0)IAXI=2
  89. C
  90. C Acquisition des connectivites FACEP
  91. C
  92. TYPE=' '
  93. CALL ACMO(MTABD,'FACEP',TYPE,MELEME)
  94. IF(TYPE.NE.'MAILLAGE')GO TO 90
  95. SEGACT MELEME
  96. C
  97. NBSOUS=LISOUS(/1)
  98. IF(NBSOUS.EQ.0)NBSOUS=1
  99. C**********************
  100. C CAS DE LA DIMENSION 2
  101. C**********************
  102. IF (IDIM.EQ.2) THEN
  103. C
  104. DO 1 L=1,NBSOUS
  105. IPT1=MELEME
  106. IF(NBSOUS.NE.1)IPT1=LISOUS(L)
  107. SEGACT IPT1
  108. NP=IPT1.NUM(/1)-1
  109. NEL=IPT1.NUM(/2)
  110. C Toutes les faces sont des SEG3 donc NP=2
  111. IF(NP.NE.2)RETURN
  112. C
  113. C On considere la face sans son centre comme un elt fini SEG2
  114. C
  115. TYPE = NOMS(IPT1.ITYPEL)//' '
  116. CALL OPTLI(IP,LIST1,TYPE,NBO)
  117. IF (IP .EQ. 0) CALL ERREUR(5)
  118. TYPE = LIST2(IP)
  119. CALL KALPBG(TYPE,'FONFORM0',IZFFM)
  120. SEGACT IZFFM*MOD
  121. IZHR=KZHR(1)
  122. SEGACT IZHR*MOD
  123. NES=GR(/1)
  124. NPG=GR(/3)
  125. C Boucle sur toutes les faces pour le paquet L
  126. DO 10 K=1,NEL
  127. NF=LECT(IPT1.NUM(NP+1,K))
  128. C REMPLISSAGE DE XYZ
  129. C ------------------
  130. C
  131. DO 12 I=1,NP
  132. J=IPT1.NUM(I,K)
  133. DO 12 N=1,IDIM
  134. XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N)
  135. 12 CONTINUE
  136. C
  137. C CALJBR calcul le jacobien du passage de l'elt de ref.
  138. C a l'elt. reel
  139. C A2J=Jacobien AIRE=detA2J
  140. CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,
  141. &NPG,IAXI,AIRE,A2J,SGEN)
  142. C
  143. NORMX=A2J(1,2,1)
  144. NORMY=A2J(2,2,1)
  145. C On verifie que n est dans le meme sens que vec(13) de FACEL
  146. C Calcul de vec(13)
  147. C FACEL est repere par IPT2
  148. J1=IPT2.NUM(1,NF)
  149. J2=IPT2.NUM(2,NF)
  150. JJ2=LECT(J2)
  151. 7 J3=IPT2.NUM(3,NF)
  152. C
  153. IF (J1.eq.J3) THEN
  154. X1=XCOOR((J1-1)*(IDIM+1)+1)
  155. Y1=XCOOR((J1-1)*(IDIM+1)+2)
  156. X2=XCOOR((J2-1)*(IDIM+1)+1)
  157. Y2=XCOOR((J2-1)*(IDIM+1)+2)
  158. SCAL=(X2-X1)*NORMX+(Y2-Y1)*NORMY
  159. SGN=SIGN(1.D0,SCAL)
  160. C
  161. ELSE
  162. X1=XCOOR((J1-1)*(IDIM+1)+1)
  163. Y1=XCOOR((J1-1)*(IDIM+1)+2)
  164. X3=XCOOR((J3-1)*(IDIM+1)+1)
  165. Y3=XCOOR((J3-1)*(IDIM+1)+2)
  166. SCAL=(X3-X1)*NORMX+(Y3-Y1)*NORMY
  167. SGN=SIGN(1.D0,SCAL)
  168. ENDIF
  169. MPOVAL.VPOCHA(JJ2,1)=SGN*NORMX
  170. MPOVAL.VPOCHA(JJ2,2)=SGN*NORMY
  171. C Calcul de la longueur
  172. C
  173. CFT=1.D0
  174. IF(IAXI.NE.0)THEN
  175. TAMP1=(XYZ(1,2)+XYZ(1,1))*0.5D0
  176. TAMP2=(XYZ(2,2)+XYZ(2,1))*0.5D0
  177. ENDIF
  178. C Axisymetrie / ox (non teste)
  179. IF (IAXI.EQ.1) CFT=2.D0*XPI*TAMP2
  180. C Axisymetrie / oy (non teste non plus)
  181. IF (IAXI.EQ.2) CFT=2.D0*XPI*TAMP1
  182. C
  183. MPOVA2.VPOCHA(JJ2,1)=AIRE*CFT
  184. C
  185. C
  186. 10 CONTINUE
  187. C
  188. SEGDES IPT1
  189. 1 CONTINUE
  190. C**********************
  191. C CAS DE LA DIMENSION 3
  192. C**********************
  193. ELSE
  194. C
  195. DO 21 L=1,NBSOUS
  196. IPT1=MELEME
  197. IF(NBSOUS.NE.1)IPT1=LISOUS(L)
  198. SEGACT IPT1
  199. NP=IPT1.NUM(/1)-1
  200. NEL=IPT1.NUM(/2)
  201. C Les elts complets ne sont pas implementes, on utilise
  202. C l'elt. face sans le centre
  203. C
  204. TYPE = NOMS(IPT1.ITYPEL)//' '
  205. CALL OPTLI(IP,LIST1,TYPE,NBO)
  206. IF (IP .EQ. 0) CALL ERREUR(5)
  207. TYPE = LIST2(IP)
  208. C Calcul des fonctions de forme sur l'elt. de reference
  209. C
  210. CALL KALPBG(TYPE,'FONFORM0',IZFFM)
  211. SEGACT IZFFM*MOD
  212. IZHR=KZHR(1)
  213. SEGACT IZHR*MOD
  214. NES=GR(/1)
  215. NPG=GR(/3)
  216. C Boucle sur les faces du paquet L
  217. DO 210 K=1,NEL
  218. NF=LECT(IPT1.NUM(NP+1,K))
  219. C
  220. C REMPLISSAGE DE XYZ
  221. C ------------------
  222. C
  223. DO 212 I=1,NP
  224. J=IPT1.NUM(I,K)
  225. DO 212 N=1,IDIM
  226. XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N)
  227. C PRINT*,XY3(N,I)
  228. 212 CONTINUE
  229. C
  230. C calcul du jacobien idem dim 2
  231. C
  232. CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,
  233. & NPG,IAXI,AIRE,A3J,SGEN)
  234. NORMX=A3J(1,3,1)
  235. NORMY=A3J(2,3,1)
  236. NORMZ=A3J(3,3,1)
  237. C
  238. C On verifie que n est dans le meme sens que vec(13) de FACEL
  239. C Calcul de vec(13)
  240. C FACEL est repere par IPT2
  241. J1=IPT2.NUM(1,NF)
  242. J2=IPT2.NUM(2,NF)
  243. JJ2=LECT(J2)
  244. J3=IPT2.NUM(3,NF)
  245. C
  246. IF (J1.eq.J3) THEN
  247. X1=XCOOR((J1-1)*(IDIM+1)+1)
  248. Y1=XCOOR((J1-1)*(IDIM+1)+2)
  249. Z1=XCOOR((J1-1)*(IDIM+1)+3)
  250. X2=XCOOR((J2-1)*(IDIM+1)+1)
  251. Y2=XCOOR((J2-1)*(IDIM+1)+2)
  252. Z2=XCOOR((J2-1)*(IDIM+1)+3)
  253. C
  254. SCAL=(X2-X1)*NORMX+(Y2-Y1)*NORMY+(Z2-Z1)*NORMZ
  255. SGN=SIGN(1.D0,SCAL)
  256. C
  257. ELSE
  258. X1=XCOOR((J1-1)*(IDIM+1)+1)
  259. Y1=XCOOR((J1-1)*(IDIM+1)+2)
  260. Z1=XCOOR((J1-1)*(IDIM+1)+3)
  261. X3=XCOOR((J3-1)*(IDIM+1)+1)
  262. Y3=XCOOR((J3-1)*(IDIM+1)+2)
  263. Z3=XCOOR((J3-1)*(IDIM+1)+3)
  264. C
  265. SCAL=(X3-X1)*NORMX+(Y3-Y1)*NORMY+(Z3-Z1)*NORMZ
  266. SGN=SIGN(1.D0,SCAL)
  267. ENDIF
  268. C
  269. MPOVAL.VPOCHA(JJ2,1)=SGN*NORMX
  270. MPOVAL.VPOCHA(JJ2,2)=SGN*NORMY
  271. MPOVAL.VPOCHA(JJ2,3)=SGN*NORMZ
  272. C
  273. C Calcul de la surface de la face
  274. C
  275. MPOVA2.VPOCHA(JJ2,1)=AIRE
  276. 210 CONTINUE
  277. C
  278. SEGDES IPT1
  279. 21 CONTINUE
  280. C
  281. ENDIF
  282. C
  283. C Orientation des normales (partie non modifiee)
  284. C
  285. SEGDES MELEMD,MELEME
  286. SEGDES MCHPOI,MPOVAL,MCHPO2,MPOVA2
  287. C
  288. TYPE=' '
  289. CALL ACMO(MTABD,'ELTFA',TYPE,MELEMT)
  290. IF(MELEMT.EQ.0)GO TO 90
  291.  
  292. TYPE=' '
  293. CALL ACMO(MTABD,'CENTRE',TYPE,MELEMC)
  294. IF(MELEMC.EQ.0)GO TO 90
  295. SEGACT MELEMC
  296.  
  297. TYPE=' '
  298. CALL ACMO(MTABD,'FACEL',TYPE,MELEMF)
  299. IF(MELEMF.EQ.0)GO TO 90
  300. SEGACT MELEMF
  301. CALL CRCHPE(MELEMT,1,MCHELM)
  302. SEGACT MCHELM
  303. NBSOUS=IMACHE(/1)
  304. C
  305. K0=0
  306. DO 41 L=1,NBSOUS
  307. MCHAML=ICHAML(L)
  308. SEGACT MCHAML
  309. MELVAL=IELVAL(1)
  310. SEGACT MELVAL*MOD
  311. IPT1=IMACHE(L)
  312. SEGACT IPT1
  313. NF=IPT1.NUM(/1)
  314. NEL=IPT1.NUM(/2)
  315. DO 42 K=1,NEL
  316. K0=K0+1
  317. NE=MELEMC.NUM(1,K0)
  318. DO 42 I=1,NF
  319. KF=IPT1.NUM(I,K)
  320. KF=LECT(KF)
  321. N1=MELEMF.NUM(1,KF)
  322. N3=MELEMF.NUM(3,KF)
  323. C
  324. IF(NE.EQ.N1)THEN
  325. C La normale est sortante
  326. A=1.D0
  327. ELSEIF(NE.EQ.N3)THEN
  328. C La normale est entrante
  329. A=-1.D0
  330. ELSE
  331. WRITE(6,*)'Problemes dans KNRF NE,N1,N3=',NE,N1,N3
  332. A=0.
  333. ENDIF
  334. C
  335. VELCHE(I,K)=A
  336. C
  337. 42 CONTINUE
  338. SEGDES MELVAL,MCHAML
  339. SEGDES IPT1
  340. 41 CONTINUE
  341. SEGSUP MLENTI
  342. SEGDES MCHELM
  343. SEGDES MELEMC,MELEMF
  344. C Normales aux faces par element
  345. C write(6,*)' faces elt Normales aire',MCHELM,MCHPOI,MCHPO2
  346. C CALL ECROBJ('MCHAML ',MCHELM)
  347. C Normales aux faces
  348. C CALL ECROBJ('CHPOINT ',MCHPOI)
  349. C aire des faces
  350. C CALL ECROBJ('CHPOINT ',MCHPO2)
  351. RETURN
  352. C
  353. 90 CONTINUE
  354. WRITE(6,*)' Interruption anormale de KNRF'
  355. RETURN
  356. END
  357.  
  358.  
  359.  
  360.  
  361.  
  362.  
  363.  
  364.  
  365.  
  366.  
  367.  
  368.  
  369.  
  370.  
  371.  
  372.  
  373.  
  374.  
  375.  
  376.  
  377.  

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