Télécharger knrf.eso

Retour à la liste

Numérotation des lignes :

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

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