Télécharger ydudw.eso

Retour à la liste

Numérotation des lignes :

ydudw
  1. C YDUDW SOURCE GOUNAND 23/08/22 21:15:04 11725
  2. SUBROUTINE YDUDW
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. C***********************************************************************
  6. C
  7. C CET OPERATEUR DISCRETISE L'OPERATEUR DE PENALISATION DIV(U)=EPS*P
  8. C EN 2D SUR LES ELEMENTS QUA4 TRI3 TRI7 et QUA9 PLAN OU AXI
  9. C EN 3D SUR LES ELEMENTS CUB8 ET PRI6
  10. C
  11. C SYNTAXE :
  12. C ---------
  13. C
  14. C DUDW(EPS) INCO UN :
  15. C
  16. C COEFFICIENTS :
  17. C --------------
  18. C
  19. C
  20. C EPS (SCAL DOMA) PARAMETRE DE PENALISATION
  21. C (SCAL ELEM)
  22. C
  23. C INCONNUES :
  24. C -----------
  25. C
  26. C UN CHAMP DE VITESSE
  27. C
  28. C
  29. C OPTIONS : POROSITE DIV(PU)=EPS*P
  30. C SOURCE DIV(U)-Q=EPS*P
  31. C
  32. C************************************************************************
  33.  
  34.  
  35. -INC PPARAM
  36. -INC CCOPTIO
  37. -INC CCGEOME
  38. -INC SMCOORD
  39. -INC SMLENTI
  40. -INC SMELEME
  41. POINTEUR MELEM1.MELEME,MELEMS.MELEME,MELEMI.MELEME,MELEMP.MELEME
  42. POINTEUR MELEMC.MELEME
  43. -INC SMCHPOI
  44. POINTEUR IZGG1.MPOVAL,SRCE.MPOVAL,IZTGG1.MPOVAL
  45.  
  46. -INC SIZFFB
  47. POINTEUR IZF1.IZFFM
  48.  
  49. -INC SMLMOTS
  50. POINTEUR LINCO.MLMOTS
  51. SEGMENT TRAV
  52. REAL*8 AF(NP,NP,NINC,NINC),AS(NP,NINC),CT(MP1,NP,NINC),PQ(MP1)
  53. ENDSEGMENT
  54. CHARACTER*8 NOMZ,TYPE,TYPC
  55. CHARACTER*(LOCOMP) NOM0,NOM1,NOM2,NOM3,NOMI
  56. DIMENSION IXV(3)
  57. PARAMETER (NTB=1)
  58. CHARACTER*8 LTAB(NTB)
  59. DIMENSION KTAB(NTB)
  60. DATA LTAB/'KIZX '/
  61. C*****************************************************************************
  62. CDUDW
  63. C write(6,*)' Operateur DUDW '
  64.  
  65. CALL LITABS(LTAB,KTAB,NTB,1,IRET)
  66. IF(IERR.NE.0)RETURN
  67. MTABX=KTAB(1)
  68. C
  69. C- Récupération de la table EQEX (pointeur MTAB1)
  70. C
  71. CALL LEKTAB(MTABX,'EQEX',MTAB1)
  72. IF(MTAB1.EQ.0)THEN
  73. C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
  74. MOTERR( 1: 8) = ' EQEX '
  75. MOTERR( 9:16) = ' EQEX '
  76. MOTERR(17:24) = ' KIZX '
  77. CALL ERREUR(786)
  78. RETURN
  79. ENDIF
  80. CALL ACME(MTAB1,'NAVISTOK',NASTOK)
  81. IF(NASTOK.EQ.0)THEN
  82. CALL ZDUDW(MTABX,MTAB1)
  83. RETURN
  84. ENDIF
  85. C
  86. C- Récupération de la table INCO (pointeur KINC)
  87. C
  88. CALL LEKTAB(MTAB1,'INCO',KINC)
  89. IF(KINC.EQ.0)THEN
  90. C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
  91. MOTERR( 1: 8) = ' INCO '
  92. MOTERR( 9:16) = ' INCO '
  93. MOTERR(17:24) = ' EQEX '
  94. CALL ERREUR(786)
  95. RETURN
  96. ENDIF
  97.  
  98. C*****************************************************************************
  99. C OPTIONS
  100. C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> SEMI
  101. C KFORM = 0 -> SI 1 -> EF 2 -> VF 3 -> EFMC
  102. C IDCEN = 0-> rien 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG
  103.  
  104. IAXI=0
  105. IF(IFOMOD.EQ.0)IAXI=2
  106. C
  107. C- Récupération de la table des options KOPT (pointeur KOPTI)
  108. C
  109. CALL LEKTAB(MTABX,'KOPT',KOPTI)
  110. IF (KOPTI.EQ.0) THEN
  111. C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
  112. MOTERR( 1: 8) = ' KOPT '
  113. MOTERR( 9:16) = ' KOPT '
  114. MOTERR(17:24) = ' KIZX '
  115. CALL ERREUR(786)
  116. RETURN
  117. ENDIF
  118.  
  119. CALL ACME(KOPTI,'KIMPL',KIMPL)
  120. CALL ACME(KOPTI,'KPOIN',KPRE)
  121. CALL ACME(KOPTI,'KFORM',KFORM)
  122.  
  123. IF(KFORM.NE.0.AND.KFORM.NE.1)THEN
  124. C Option %m1:8 incompatible avec les données
  125. MOTERR( 1: 8) = 'EF/EFM1 '
  126. CALL ERREUR(803)
  127. RETURN
  128. ENDIF
  129.  
  130. C write(6,*)' Apres les options '
  131. C*****************************************************************************
  132. C
  133. C- Récupération de la table DOMAINE associée au domaine local
  134. C
  135. CALL ACMM(MTABX,'NOMZONE',NOMZ)
  136. CALL LEKTAB(MTABX,'DOMZ',MTABZ)
  137. IF(MTABZ.EQ.0)THEN
  138. C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
  139. MOTERR( 1: 8) = ' DOMZ '
  140. MOTERR( 9:16) = ' DOMZ '
  141. MOTERR(17:24) = ' KIZX '
  142. CALL ERREUR(786)
  143. RETURN
  144. ENDIF
  145.  
  146. CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
  147. CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
  148. CALL LEKTAB(MTABZ,'MACRO1',MACRO1)
  149. MELEMI=MACRO1
  150. CALL LEKTAB(MTABZ,'QUADRATI',MQUAD)
  151. IF (IERR.NE.0) RETURN
  152.  
  153. IF(MQUAD.EQ.0.AND.MACRO1.EQ.0.AND.KPRE.NE.2)THEN
  154. WRITE(6,*)'Operateur DUDW '
  155. WRITE(6,*)'Type d''éléments non prévu'
  156. RETURN
  157. ENDIF
  158.  
  159. IF(KPRE.EQ.2)THEN
  160. CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
  161. MELEMP=MELEMC
  162. MELEMI=MELEME
  163. ELSEIF(KPRE.EQ.3)THEN
  164. CALL LEKTAB(MTABZ,'CENTREP0',MELEMC)
  165. MELEMP=MELEMC
  166. ELSEIF(KPRE.EQ.4)THEN
  167. CALL LEKTAB(MTABZ,'CENTREP1',MELEMC)
  168. CALL LEKTAB(MTABZ,'ELTP1NC ',MELEMP)
  169. ENDIF
  170.  
  171. SEGACT MELEME
  172. SEGACT MELEMP
  173. C*************************************************************************
  174. C VERIFICATIONS SUR LES INCONNUES
  175. C write(6,*)' Verification des inconnues '
  176.  
  177. TYPE='LISTMOTS'
  178. CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO)
  179. SEGACT LINCO
  180.  
  181. NBINC=LINCO.MOTS(/2)
  182. IF(NBINC.NE.1)THEN
  183. C Indice %m1:8 : contient plus de %i1 %m9:16
  184. MOTERR( 1:8) = 'LISTINCO'
  185. INTERR(1) = 1
  186. MOTERR(9:16) = ' MOTS '
  187. CALL ERREUR(799)
  188. RETURN
  189. ENDIF
  190.  
  191. NOMI=LINCO.MOTS(1)
  192.  
  193. TYPE=' '
  194. CALL ACMO(KINC,NOMI,TYPE,MCHPOI)
  195. IF(TYPE.NE.'CHPOINT ')THEN
  196. C Indice %m1:8 : ne contient pas un objet de type %m9:16
  197. MOTERR( 1: 8) = 'INC '//NOMI
  198. MOTERR( 9:16) = 'CHPOINT '
  199. CALL ERREUR(800)
  200. RETURN
  201. ELSE
  202. CALL LICHT(MCHPOI,IZTU1,TYPC,MELEM1)
  203. ENDIF
  204.  
  205. C*************************************************************************
  206. C Le domaine de definition est donne par le SPG de la premiere inconnue
  207. C Les inconnues suivantes devront posseder ce meme pointeur
  208. C On verifie que les points de la zone sont tous inclus dans ce SPG
  209.  
  210. CALL KRIPAD(MELEM1,MLENTI)
  211.  
  212. CALL VERPAD(MLENTI,MELEME,IRET)
  213. IF(IRET.NE.0)THEN
  214. C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
  215. MOTERR(1: 8) = 'INC '//NOMI
  216. MOTERR(9:16) = 'CHPOINT '
  217. CALL ERREUR(788)
  218. RETURN
  219. ENDIF
  220.  
  221. NINC=IDIM
  222.  
  223. MLENT1=MLENTI
  224.  
  225. C*****************************************************************************
  226. C Lecture du ou des coefficients
  227. C write(6,*)' Lecture des coefficients'
  228.  
  229. CALL ACME(MTABX,'IARG',IARG)
  230.  
  231. C? IXV(1)=MELEMC
  232. IXV(1)=0
  233. IXV(2)=1
  234. IXV(3)=0
  235. CALL LEKCOF('Opérateur DUDW :',
  236. & MTABX,KINC,1,IXV,IZTG1,IZTGG1,NPT1,NC1,IK1,IRET)
  237. IF(IRET.EQ.0)RETURN
  238. KITT=2
  239.  
  240. INDGS=0
  241. SRCE=IZTGG1
  242. MPOVA1 =IZTGG1
  243. MPOVA3 =IZTGG1
  244.  
  245. IF(IARG.EQ.2)THEN
  246. IXV(1)=MELEMC
  247. IXV(2)=1
  248. IXV(3)=0
  249. CALL LEKCOF('Opérateur DUDW :',
  250. & MTABX,KINC,2,IXV,IZS,SRCE,NPTS,NCS,IKS,IRET)
  251. IF(IRET.EQ.0)RETURN
  252. INDGS=1
  253. CALL KRIPAD(MELEMC,MLENT1)
  254. ENDIF
  255.  
  256. C*************************************************************************
  257. IF(KFORM.EQ.0)THEN
  258. C CAS FORMULATION EF SI (GRESHO)
  259.  
  260. IF(KIMPL.NE.0)THEN
  261. GO TO 90
  262. ENDIF
  263.  
  264. WRITE(6,*)' Operateur DUDW '
  265. WRITE(6,*)' Cas Formulation EF SI '
  266. WRITE(6,*)' Cas invalide '
  267. GO TO 90
  268.  
  269. C***********************************************************************
  270. ELSE IF(KFORM.EQ.1)THEN
  271. C CAS FORMULATION EF
  272.  
  273. NUTOEL=0
  274.  
  275. SEGACT MELEMI
  276. NBSOUS=MELEMI.LISOUS(/1)
  277. IF(NBSOUS.EQ.0)NBSOUS=1
  278.  
  279. NRIGE=7
  280. NKID =9
  281. NKMT =7
  282. NMATRI=1
  283. SEGINI MATRIK
  284. IRIGEL(1,1)=MELEMI
  285. IRIGEL(2,1)=MELEMI
  286. NBOP=0
  287. NBME=NINC*NINC
  288. C NBME=2
  289. NBELC=0
  290. SEGINI IMATRI
  291. IRIGEL(4,1)=IMATRI
  292. KSPGP=MELEMS
  293. KSPGD=MELEMS
  294.  
  295. IF(NBME.EQ.4)THEN
  296. NOM1='1'//NOMI(1:7)
  297. NOM2='2'//NOMI(1:7)
  298. LISPRI(1)=NOM1
  299. LISDUA(1)=NOM1
  300. LISPRI(2)=NOM2
  301. LISDUA(2)=NOM1
  302. LISPRI(3)=NOM1
  303. LISDUA(3)=NOM2
  304. LISPRI(4)=NOM2
  305. LISDUA(4)=NOM2
  306. ELSEIF(NBME.EQ.9)THEN
  307. NOM1='1'//NOMI(1:7)
  308. NOM2='2'//NOMI(1:7)
  309. NOM3='3'//NOMI(1:7)
  310. LISPRI(1)=NOM1
  311. LISDUA(1)=NOM1
  312. LISPRI(2)=NOM2
  313. LISDUA(2)=NOM1
  314. LISPRI(3)=NOM3
  315. LISDUA(3)=NOM1
  316.  
  317. LISPRI(4)=NOM1
  318. LISDUA(4)=NOM2
  319. LISPRI(5)=NOM2
  320. LISDUA(5)=NOM2
  321. LISPRI(6)=NOM3
  322. LISDUA(6)=NOM2
  323.  
  324. LISPRI(7)=NOM1
  325. LISDUA(7)=NOM3
  326. LISPRI(8)=NOM2
  327. LISDUA(8)=NOM3
  328. LISPRI(9)=NOM3
  329. LISDUA(9)=NOM3
  330.  
  331. ELSEIF(NBME.EQ.2)THEN
  332.  
  333. LISPRI(1)='1'//NOMI(1:7)
  334. LISDUA(1)='1'//NOMI(1:7)
  335. LISPRI(2)='2'//NOMI(1:7)
  336. LISDUA(2)='2'//NOMI(1:7)
  337.  
  338. ELSE
  339. WRITE(6,*)' Operateur DUDW '
  340. WRITE(6,*)' Cas invalide'
  341. GO TO 90
  342. ENDIF
  343.  
  344.  
  345. IF(INDGS.NE.0)THEN
  346. NAT=2
  347. NSOUPO=1
  348. SEGACT MELEMS
  349. N=MELEMS.NUM(/2)
  350. NC=NINC
  351. SEGINI MCHPO1,MSOUP1,MPOVA1
  352. MCHPO1.IFOPOI=IFOUR
  353. MCHPO1.MOCHDE=TITREE
  354. MCHPO1.MTYPOI='SMBR'
  355. MCHPO1.JATTRI(1)=2
  356. MCHPO1.IPCHP(1)=MSOUP1
  357. MSOUP1.NOCOMP(1)=NOM1
  358. MSOUP1.NOCOMP(2)=NOM2
  359. IF(NINC.EQ.3)MSOUP1.NOCOMP(3)=NOM3
  360. MSOUP1.IGEOC=MELEMS
  361. MSOUP1.IPOVAL=MPOVA1
  362.  
  363. SEGACT MELEMC
  364. N=MELEMC.NUM(/2)
  365. NC=1
  366. SEGINI MCHPO3,MSOUP3,MPOVA3
  367. MCHPO3.IFOPOI=IFOUR
  368. MCHPO3.MOCHDE=TITREE
  369. MCHPO3.MTYPOI='SMBR'
  370. MCHPO3.JATTRI(1)=2
  371. MCHPO3.IPCHP(1)=MSOUP3
  372. MSOUP3.NOCOMP(1)='SCAL'
  373. MSOUP3.IGEOC=MELEMC
  374. MSOUP3.IPOVAL=MPOVA3
  375.  
  376. ENDIF
  377.  
  378. DO 101 L=1,NBSOUS
  379. IPT1=MELEMI
  380. IPT2=MELEMI
  381. IF(NBSOUS.NE.1)IPT1=MELEMI.LISOUS(L)
  382. SEGACT IPT1
  383. IF(INDGS.NE.0)THEN
  384. IPT2=MELEMP
  385. IF(NBSOUS.NE.1)IPT2=MELEMP.LISOUS(L)
  386. SEGACT IPT2
  387. ENDIF
  388.  
  389. IF(MQUAD.NE.0)THEN
  390. IF(KPRE.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//'PRP0'
  391. IF(KPRE.EQ.3)NOM0=NOMS(IPT1.ITYPEL)//'PRP0'
  392. IF(KPRE.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'PRP1'
  393. ELSEIF(MACRO1.NE.0)THEN
  394. IF(KPRE.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//' '
  395. IF(KPRE.EQ.3)NOM0=NOMS(IPT1.ITYPEL)//'MCP0'
  396. IF(KPRE.EQ.4)NOM0=NOMS(IPT1.ITYPEL)//'MCP1'
  397. ELSE
  398. IF(KPRE.EQ.2)NOM0=NOMS(IPT1.ITYPEL)//' '
  399. ENDIF
  400.  
  401. CALL KALPBG(NOM0,'FONFORM ',IZFFM)
  402.  
  403.  
  404. SEGACT IZFFM*MOD
  405. IZHR=KZHR(1)
  406. SEGACT IZHR*MOD
  407. NES=GR(/1)
  408. NPG=GR(/3)
  409. IZF1=KTP(1)
  410. SEGACT IZF1*MOD
  411. MP1=IZF1.FN(/1)
  412. NP = IPT1.NUM(/1)
  413.  
  414. MP = NP
  415. NBEL=IPT1.NUM(/2)
  416. SEGINI IPM1,IPM2,IPM3,IPM4
  417. LIZAFM(L,1)=IPM1
  418. C LIZAFM(L,2)=IPM4
  419. LIZAFM(L,2)=IPM2
  420. LIZAFM(L,3)=IPM3
  421. LIZAFM(L,4)=IPM4
  422. IPM5=IPM1
  423. IPM6=IPM1
  424. IPM7=IPM1
  425. IPM8=IPM1
  426. IPM9=IPM1
  427. IF(NBME.EQ.9)THEN
  428. SEGINI IPM5,IPM6,IPM7,IPM8,IPM9
  429. LIZAFM(L,5)=IPM5
  430. LIZAFM(L,6)=IPM6
  431. LIZAFM(L,7)=IPM7
  432. LIZAFM(L,8)=IPM8
  433. LIZAFM(L,9)=IPM9
  434. ENDIF
  435.  
  436. C Pour l'instant les sources et puits de masse ne sont pas actifs
  437. INDG2=0
  438. KI2=0
  439. KJ2=0
  440.  
  441. SEGINI TRAV
  442.  
  443. NPT=MPOVA1.VPOCHA(/1)
  444. SEGACT,MCOORD
  445. CALL XDUDW(FN,IZF1.FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ,
  446. & NES,IDIM,NP,MP1,NPG,IAXI,NINC,
  447. & IZTGG1.VPOCHA,IK1,SRCE.VPOCHA,INDGS,IKS,
  448. & IPT1.NUM,NBEL,NUTOEL,XCOOR,AF,AS,CT,PQ,
  449. & IPM1.AM,IPM2.AM,IPM3.AM,IPM4.AM,IPM5.AM,IPM6.AM,IPM7.AM,
  450. & IPM8.AM,IPM9.AM,MPOVA1.VPOCHA,NPT,LECT,IPT2.NUM,MLENT1.LECT,
  451. & MPOVA3.VPOCHA)
  452. SEGDES,MCOORD
  453.  
  454. SEGSUP TRAV
  455. NUTOEL=NUTOEL+NBEL
  456. 101 CONTINUE
  457. SEGDES IMATRI,MATRIK
  458.  
  459. CALL ECROBJ('MATRIK',MATRIK)
  460. C? CALL ECMO(MTABX,'MATELM','MATRIK',MATRIK)
  461.  
  462. NAT=2
  463. NSOUPO=0
  464. SEGINI MCHPOI
  465. JATTRI(1)=2
  466. CALL ECROBJ('CHPOINT',MCHPOI)
  467.  
  468.  
  469.  
  470. IF(INDGS.NE.0)THEN
  471.  
  472. TYPE=' '
  473. CALL ACMO(MTAB1,'SMBR',TYPE,MCHPO2)
  474. IF(TYPE.NE.'CHPOINT')THEN
  475. CALL ECMO(MTAB1,'SMBR','CHPOINT',MCHPO1)
  476. ELSE
  477. CALL ECROBJ('CHPOINT',MCHPO2)
  478. CALL ECROBJ('CHPOINT',MCHPO1)
  479. CALL PRFUSE
  480. CALL LIROBJ('CHPOINT',MCHPOI,1,IRET)
  481. CALL ECMO(MTAB1,'SMBR','CHPOINT',MCHPOI)
  482. ENDIF
  483.  
  484. TYPE=' '
  485. C CALL ACMO(MTAB1,'DUDWSRCE',TYPE,MCHPO4)
  486. C IF(TYPE.NE.'CHPOINT')THEN
  487. CALL ECMO(MTAB1,'DUDWSRCE','CHPOINT',MCHPO3)
  488. C ELSE
  489. C CALL ECROBJ('CHPOINT',MCHPO3)
  490. C CALL ECROBJ('CHPOINT',MCHPO4)
  491. C CALL PRFUSE
  492. C CALL LIROBJ('CHPOINT',MCHPOI,1,IRET)
  493. C CALL ECMO(MTAB1,'DUDWSRCE','CHPOINT',MCHPOI)
  494. C ENDIF
  495.  
  496. ENDIF
  497.  
  498. SEGSUP MLENTI
  499. IF(INDGS.NE.0)SEGSUP MLENT1
  500.  
  501. RETURN
  502. C*************************************************************************
  503. ELSE IF(KFORM.EQ.2)THEN
  504. C CAS FORMULATION VF
  505.  
  506. RETURN
  507. ENDIF
  508. C*************************************************************************
  509. 1002 FORMAT(10(1X,1PE11.4))
  510. 90 CONTINUE
  511. * WRITE(6,*)' Interuption anormale de DUDW'
  512. * 34 2
  513. * Operateur incompatible avec les options utilisees
  514. CALL ERREUR(34)
  515. RETURN
  516. END
  517.  
  518.  

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