Télécharger knol.eso

Retour à la liste

Numérotation des lignes :

  1. C KNOL SOURCE MAGN 17/02/24 21:15:16 9323
  2. SUBROUTINE KNOL
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. C*************************************************************************
  6. C Operateur KNOL
  7. C
  8. C Objet Transforme un CHAMPOINT SOMMET en un CHAMPOINT CENTRE
  9. C
  10. C SYNTAXE : CHPC = KNOL TABDOM CHPS ;
  11. C
  12. C TABDOM : Table DOMAINE contenant les supports geometriques de CHPS
  13. C CHPS : CHAMPOINT SOMMET
  14. C CHPC : CHAMPOINT CENTRE
  15. C
  16. C
  17. C*************************************************************************
  18. -INC CCOPTIO
  19. -INC CCREEL
  20. -INC CCGEOME
  21. -INC SMCOORD
  22. -INC SMELEME
  23. POINTEUR MELEMS.MELEME,MELEMD.MELEME,SPGD.MELEME
  24. -INC SMMODEL
  25. -INC SMCHPOI
  26. POINTEUR IZB.MCHPOI,IZBB.MPOVAL
  27. -INC SMLENTI
  28. -INC SIZFFB
  29. POINTEUR IZF1.IZFFM,IZH2.IZHR,IZW.IZFFM,IZWH.IZHR
  30. SEGMENT SAJT
  31. REAL*8 AJT(IDIM,IDIM,NPG),RF1(NP,MP,IDIM),SM1(NP,IDIM)
  32. c REAL*8 TN1(NP,IDIM),TN2(NP,IDIM)
  33. ENDSEGMENT
  34.  
  35. PARAMETER (NTO=4,NBMO=4)
  36. DIMENSION ITABO(NTO)
  37. CHARACTER*4 NOMD4
  38. CHARACTER*8 TYPE,TYPC,LISMO(NBMO),TYPSPG,MTERR,NOM0
  39. DATA LISMO/'CENTRE ',' ', 'CENTREP1','MSOMMET'/
  40. C***
  41.  
  42. TYPSPG='CENTRE '
  43. CALL INITI(ITABO,NTO,0)
  44. 4 CONTINUE
  45.  
  46. CALL QUETYP(TYPE,0,IRET)
  47. c write(6,*)' Iret de Quetyp= ',iret,' TYPSPG=',typspg
  48.  
  49. IF(IRET.EQ.0)THEN
  50. IF(TYPSPG.EQ.'MSOMMET '.OR.TYPSPG.EQ.'CENTRE '.OR.
  51. & TYPSPG.EQ.'CENTREP1'.OR.TYPSPG.EQ.' ')THEN
  52. IF(ITABO(1).EQ.1.AND.ITABO(2).EQ.1)GO TO 52
  53.  
  54. IF(ITABO(1).NE.1)THEN
  55. C% Il faut spécifier un objet de type %m1:8 et de sous type %m9:16
  56. MOTERR(1: 8) = 'CHPOINT '
  57. MOTERR(9:16) = 'DIFFUS '
  58. CALL ERREUR(79)
  59. ENDIF
  60.  
  61. IF(ITABO(2).NE.1)THEN
  62. C% Il faut spécifier un objet de type %m1:8 et de sous type %m9:16
  63. MOTERR(1: 8) = 'MMODEL '
  64. MOTERR(9:16) = ' '
  65. CALL ERREUR(79)
  66. ENDIF
  67.  
  68. ENDIF
  69. RETURN
  70. ENDIF
  71. *
  72. * Lecture du CHPOIN
  73. *
  74. IF(TYPE.EQ.'CHPOINT')THEN
  75. IF(ITABO(1).NE.0)THEN
  76. C% On a déja lu un objet de type %m1:8
  77. MOTERR(1: 8) = 'CHPOINT '
  78. CALL ERREUR(966)
  79. RETURN
  80. ENDIF
  81. ITABO(1)=1
  82.  
  83. CALL LIROBJ('CHPOINT',IZB,1,IRETOU)
  84. SEGACT IZB
  85. IF(IZB.IPCHP(/1).NE.1)THEN
  86. C% Erreur dans le partitionnement
  87. CALL ERREUR(920)
  88. SEGDES IZB
  89. RETURN
  90. ENDIF
  91. GO TO 4
  92. ENDIF
  93. *
  94. * Lecture de l'objet modele 'Navier-Stokes'
  95. *
  96. C***
  97. IF(TYPE.EQ.'MMODEL ')THEN
  98. IF(ITABO(2).NE.0)THEN
  99. C% On a déja lu un objet de type %m1:8
  100. MOTERR(1: 8) = 'MMODEL '
  101. MOTERR(9:16) = ' '
  102. CALL ERREUR(966)
  103. RETURN
  104. ENDIF
  105. ITABO(2)=1
  106.  
  107. CALL LIROBJ('MMODEL ',MMODEL,1,IRET)
  108. SEGACT MMODEL
  109. N1=KMODEL(/1)
  110. DO 41 L=1,N1
  111. IMODEL=KMODEL(L)
  112. SEGACT IMODEL
  113. IF(FORMOD(1).NE.'NAVIER_STOKES')THEN
  114. IF(FORMOD(1).NE.'DARCY')THEN
  115. C% On veut un modèle de type %m1:16 .
  116. MOTERR( 1:16) = 'NAVIER_STOKES '
  117. CALL ERREUR(719)
  118. RETURN
  119. ENDIF
  120. ENDIF
  121. SEGDES IMODEL
  122. 41 CONTINUE
  123.  
  124. CALL LEKMOD(MMODEL,MTABD,INEFMD)
  125. C /S INEFMD : Type formulation INEFMD=1 LINE,=2 MACRO,=3 QUADRATIQUE
  126. C INEFMD=4 LINB
  127. C KPOIN = 0->SOMMET 1-> FACE 2-> CENTRE 3-> CENTREP0 4-> CENTREP1 5-> MSOMMET
  128. GO TO 4
  129. ENDIF
  130. *
  131. * Lecture de l'objet table DOMAINE
  132. *
  133. C***
  134. IF(TYPE.EQ.'TABLE ')THEN
  135. IF(ITABO(4).NE.0)THEN
  136. C% On a déja lu un objet de type %m1:8
  137. MOTERR(1: 8) = 'TABLE '
  138. MOTERR(9:16) = ' '
  139. CALL ERREUR(966)
  140. RETURN
  141. ENDIF
  142. ITABO(2)=1
  143.  
  144. CALL LIROBJ('TABLE ',MTABD,1,IRET)
  145. GO TO 4
  146. ENDIF
  147. *
  148. * Lecture d'un mot
  149. *
  150. IF(TYPE.EQ.'MOT ')THEN
  151. CALL LIRMOT(LISMO,NBMO,IP,1)
  152. TYPSPG=LISMO(IP)
  153. GO TO 4
  154. ENDIF
  155.  
  156.  
  157.  
  158. 52 CONTINUE
  159. C
  160. C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  161. C------------------ Traitement du cas MSOMMET ---------------------------
  162.  
  163. IF(TYPSPG.EQ.'MSOMMET')THEN
  164. CALL LICHTL(IZB,IZBB,TYPC,IGEOMB)
  165. NC=IZBB.VPOCHA(/2)
  166. TYPE=' '
  167. CALL ACMO(MTABD,'SOMMET',TYPE,MELEMS)
  168. IF(TYPE.NE.'MAILLAGE')GO TO 90
  169.  
  170. C
  171. C On verifie Le support du CHP1 (SOMMET)
  172. C
  173. CALL KRIPAD(IGEOMB,MLENTI)
  174. CALL VERPAD(MLENTI,MELEMS,IRET)
  175.  
  176. IF(IRET.NE.0)THEN
  177. C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
  178. MOTERR(1: 8) = 'SOMMET'
  179. MOTERR(9:16) = 'CHPOINT '
  180. SEGSUP MLENTI
  181. RETURN
  182. ENDIF
  183. SEGSUP MLENTI
  184. SEGDES MELEMS
  185.  
  186. TYPE=' '
  187. CALL ACMO(MTABD,'MSOMMET',TYPE,MELEME)
  188. IF(TYPE.NE.'MAILLAGE')GO TO 90
  189.  
  190. CALL ECROBJ('MAILLAGE',MELEME)
  191. CALL ECROBJ('CHPOINT',IZB)
  192. CALL REDU
  193.  
  194. RETURN
  195. ENDIF
  196. C-------------- FIN Traitement du cas MSOMMET ---------------------------
  197. C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  198.  
  199. C
  200. C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  201. C------------------ Traitement du cas CENTRE ----------------------------
  202.  
  203. IF(TYPSPG.EQ.'CENTRE')THEN
  204. CALL LICHTL(IZB,IZBB,TYPC,IGEOMB)
  205. NC=IZBB.VPOCHA(/2)
  206. TYPE=' '
  207. CALL ACMO(MTABD,'SOMMET',TYPE,MELEMS)
  208. IF(TYPE.NE.'MAILLAGE')GO TO 90
  209.  
  210.  
  211.  
  212. CALL KRIPAD(IGEOMB,MLENTI)
  213. CALL VERPAD(MLENTI,MELEMS,IRET)
  214.  
  215. IF(IRET.NE.0)THEN
  216. C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
  217. MOTERR(1: 8) = 'SOMMET'
  218. MOTERR(9:16) = 'CHPOINT '
  219. SEGSUP MLENTI
  220. RETURN
  221. ENDIF
  222.  
  223. TYPE=' '
  224. CALL ACMO(MTABD,'CENTRE',TYPE,MELEMC)
  225. IF(TYPE.NE.'MAILLAGE')GO TO 90
  226.  
  227. TYPE=' '
  228. CALL ACMO(MTABD,'MAILLAGE',TYPE,MELEME)
  229. IF(TYPE.NE.'MAILLAGE')GO TO 90
  230.  
  231.  
  232. TYPE='CENTRE'
  233. c write(6,*)' Apparemment on traite le cas centre !!!'
  234. CALL CRCHPT(TYPE,MELEMC,NC,MCHPOI)
  235. CALL LICHTM(MCHPOI,MPOVAL,TYPC,IGEOM)
  236.  
  237. SEGACT MELEME
  238.  
  239. NBSOUS=LISOUS(/1)
  240. IF(NBSOUS.EQ.0)NBSOUS=1
  241.  
  242. KK=0
  243. DO 1 L=1,NBSOUS
  244. IPT1=MELEME
  245. IF(NBSOUS.NE.1)IPT1=LISOUS(L)
  246. SEGACT IPT1
  247. NP=IPT1.NUM(/1)
  248. NEL=IPT1.NUM(/2)
  249. DO 10 K=1,NEL
  250. KK=KK+1
  251. DO 13 N=1,NC
  252. UU=0.D0
  253. DO 12 I=1,NP
  254. IU=LECT(IPT1.NUM(I,K))
  255. UU=UU+IZBB.VPOCHA(IU,N)
  256. 12 CONTINUE
  257. VPOCHA(KK,N)=UU/FLOAT(NP)
  258. 13 CONTINUE
  259.  
  260. 10 CONTINUE
  261. SEGDES IPT1
  262. 1 CONTINUE
  263. SEGDES IZB,IZBB
  264. SEGDES MELEME,MPOVAL,MCHPOI
  265. C
  266. SEGSUP MLENTI
  267. CALL ECROBJ('CHPOINT ',MCHPOI)
  268.  
  269. RETURN
  270. ENDIF
  271. C-------------- FIN Traitement du cas CENTRE ----------------------------
  272. C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  273.  
  274. C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  275. C------------------ Traitement du cas CENTREP1 --------------------------
  276.  
  277. IF(TYPSPG.EQ.'CENTREP1')THEN
  278.  
  279. MTERR='EF CTRP1'
  280. IF(INEFMD.EQ.2)NOMD4='MCP1'
  281. IF(INEFMD.EQ.3)NOMD4='PRP1'
  282. IF(INEFMD.NE.2.AND.INEFMD.NE.3)THEN
  283. C Option %m1:8 incompatible avec les données
  284. MOTERR( 1: 8) = MTERR
  285. CALL ERREUR(803)
  286. RETURN
  287. ENDIF
  288.  
  289. CALL LEKTAB(MTABD,'CENTREP1',SPGD)
  290. CALL LEKTAB(MTABD,'ELTP1NC ',MELEMD)
  291.  
  292.  
  293. CALL LICHTL(IZB,IZBB,TYPC,IGEOMB)
  294. NC=IZBB.VPOCHA(/2)
  295. TYPE=' '
  296. CALL ACMO(MTABD,'SOMMET',TYPE,MELEMS)
  297. IF(TYPE.NE.'MAILLAGE')GO TO 90
  298.  
  299.  
  300. CALL KRIPAD(IGEOMB,MLENT1)
  301. CALL VERPAD(MLENT1,MELEMS,IRET)
  302.  
  303. IF(IRET.NE.0)THEN
  304. C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
  305. MOTERR(1: 8) = 'SOMMET'
  306. MOTERR(9:16) = 'CHPOINT '
  307. SEGSUP MLENTI
  308. RETURN
  309. ENDIF
  310.  
  311. TYPE=' '
  312. CALL ACMO(MTABD,'MAILLAGE',TYPE,MELEME)
  313. IF(TYPE.NE.'MAILLAGE')GO TO 90
  314.  
  315. IF(INEFMD.EQ.2)CALL LEKTAB(MTABD,'MACRO1',MELEME)
  316. IF (IERR.NE.0) RETURN
  317.  
  318. TYPE='CENTREP1'
  319. c write(6,*)' Apparemment on traite le cas centrep1 !!!'
  320. CALL CRCHPT(TYPE,SPGD,NC,MCHPOI)
  321. CALL LICHTM(MCHPOI,MPOVAL,TYPC,IGEOM)
  322. c..........................................................................
  323.  
  324.  
  325. c IK3=0
  326. IAXI=0
  327. IF(IFOMOD.EQ.0)IAXI=2
  328. DEUPI=1.D0
  329. IF(IAXI.NE.0)DEUPI=2.D0*XPI
  330.  
  331. NC=MPOVAL.VPOCHA(/2)
  332. CALL KRIPAD(SPGD,MLENTI)
  333.  
  334. SEGACT MELEME
  335.  
  336. NKD=0
  337. DO 101 L=1,MAX(1,LISOUS(/1))
  338. SEGACT MELEMD
  339. IPT1=MELEME
  340. IPT2=MELEMD
  341. IF(LISOUS(/1).NE.0)IPT1=LISOUS(L)
  342. SEGACT IPT1
  343. IF(MELEMD.LISOUS(/1).NE.0)THEN
  344. IPT2=MELEMD.LISOUS(L)
  345. NKD=0
  346. ENDIF
  347. SEGACT IPT2
  348. MP=IPT2.NUM(/1)
  349.  
  350. NOM0 = NOMS(IPT1.ITYPEL)//NOMD4
  351. c write(6,*)' KNOL 1er KALPBG NOM0=',NOM0,IPT1
  352. CALL KALPBG(NOM0,'FONFORM ',IZFFM)
  353. IF(IZFFM.EQ.0)RETURN
  354. SEGACT IZFFM*MOD
  355. IZHR=KZHR(1)
  356. SEGACT IZHR*MOD
  357. IZF1 = KTP(1)
  358. IZH2 = KZHR(2)
  359. IZW = IZF1
  360. SEGACT IZW*MOD
  361. IF(MP.NE.IZW.FN(/1))THEN
  362. write(6,*)' Gross problem ds KNOL '
  363. write(6,*)' NOM0=',NOM0 ,' NOMD4=',NOMD4
  364. write(6,*)' MP=',MP,' IZW.FN(/1)='
  365. & ,IZW.FN(/1)
  366. return
  367. ENDIF
  368.  
  369. NES=GR(/1)
  370. NPG=GR(/3)
  371. NP = IPT1.NUM(/1)
  372. NBEL=IPT1.NUM(/2)
  373. SEGINI SAJT
  374. c write(6,*)' AVANT 108 NC=',NC,' NBEL=',NBEL,MP,NP,NC
  375. c write(6,*)' AVANT 108 IK4=',IK4
  376. DO 108 KE=1,NBEL
  377.  
  378. NKD=NKD+1
  379.  
  380. DO 109 I=1,NP
  381. J=IPT1.NUM(I,KE)
  382. DO 109 N=1,IDIM
  383. XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N)
  384. 109 CONTINUE
  385.  
  386. CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,
  387. * IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
  388.  
  389. CALL INITD(SM1,(MP*IDIM),0.D0)
  390.  
  391. C=======================================================================
  392.  
  393. C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  394. C...... Source
  395. DO 710 I=1,MP
  396. U2=0.D0
  397. U4=0.D0
  398. DO 717 N=1,NC
  399. DO 715 LG=1,NPG
  400. WT=IZW.FN(I,LG)
  401. U4=U4+WT*PGSQ(LG)*DEUPI*RPG(LG)
  402. UJ=0.D0
  403. DO 714 J=1,NP
  404. JU=MLENT1.LECT(IPT1.NUM(J,KE))
  405. UJ=UJ+FN(J,LG)*IZBB.VPOCHA(JU,N)
  406. 714 CONTINUE
  407. U2=U2+UJ*WT*PGSQ(LG)*DEUPI*RPG(LG)
  408. 715 CONTINUE
  409.  
  410. SM1(I,N)=SM1(I,N)+(U2/U4)
  411. 717 CONTINUE
  412. 710 CONTINUE
  413. c write(6,*)' SM1 *****'
  414. c do 711 n=1,nc
  415. c write(6,1002)(sm1(i,n),i=1,mp)
  416. c711 continue
  417. C...... Source Fin
  418. C=======================================================================
  419.  
  420.  
  421. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  422. C ...... Chargement Second membre
  423. DO 910 I=1,MP
  424. I1=LECT(IPT2.NUM(I,NKD))
  425. DO 910 N=1,NC
  426. VPOCHA(I1,N)=VPOCHA(I1,N)+SM1(I,N)
  427. 910 CONTINUE
  428. C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  429.  
  430. 108 CONTINUE
  431.  
  432.  
  433. SEGDES IPT1,IPT2
  434.  
  435. SEGSUP IZFFM,IZHR,IZF1,IZH2
  436. SEGSUP SAJT
  437.  
  438. 101 CONTINUE
  439. SEGSUP MLENTI
  440.  
  441. c..........................................................................
  442. SEGSUP MLENTI,MLENT1
  443. CALL ECROBJ('CHPOINT ',MCHPOI)
  444.  
  445. RETURN
  446. ENDIF
  447. C-------------- FIN Traitement du cas CENTREP1 --------------------------
  448. C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  449.  
  450.  
  451. 90 CONTINUE
  452. WRITE(6,*)'Interruption anormale de KLNO '
  453. RETURN
  454.  
  455. 1001 FORMAT(20(1X,I5))
  456. 1002 FORMAT(10(1X,1PE11.4))
  457. 1008 FORMAT(10(1X,A8))
  458. END
  459.  
  460.  
  461.  
  462.  
  463.  
  464.  
  465.  
  466.  
  467.  
  468.  
  469.  
  470.  
  471.  
  472.  
  473.  
  474.  
  475.  
  476.  

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