Télécharger divs.eso

Retour à la liste

Numérotation des lignes :

divs
  1. C DIVS SOURCE CB215821 20/11/25 13:25:40 10792
  2. SUBROUTINE DIVS(IKAS,MTABX)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. C************************************************************************
  6. C Operateur KMAC
  7. C
  8. C OBJET : Cree un objet de type MATRIK
  9. C
  10. C
  11. C IKAS=1 KMAB calcul de B (DIV U)
  12. C IKAS=2 KMBT calcul de Bt uniquement (GRAD P)
  13. C IKAS=3 KBBT calcul de B assemblage pour B et Bt
  14. C
  15. C***********************************************************************
  16.  
  17. -INC PPARAM
  18. -INC CCOPTIO
  19. -INC CCGEOME
  20. -INC CCREEL
  21. -INC SMCOORD
  22. -INC SIZFFB
  23. -INC SMCHPOI
  24. POINTEUR IZCH2.MCHPOI,IZCCH2.MPOVAL
  25. POINTEUR IZDV.MCHPOI,IZDDV.MPOVAL,IZTU1.MPOVAL,TETAN.MPOVAL
  26. POINTEUR IZTG1.MCHPOI,IZTGG1.MPOVAL,IZBETA.MPOVAL
  27. POINTEUR IZTGG.MCHPOI,IZANU.MPOVAL
  28.  
  29. -INC SMLENTI
  30. POINTEUR IZIPAD.MLENTI,MLENTI1.MLENTI,MLENTI2.MLENTI,IPADU.MLENTI
  31. -INC SMLMOTS
  32. POINTEUR LINCO.MLMOTS
  33. -INC SMELEME
  34. POINTEUR MELEMZ.MELEME,MELEMB.MELEME,MELSTB.MELEME
  35. POINTEUR MELEM1.MELEME,MELES1.MELEME,MCTREI.MELEME
  36. POINTEUR IGEOM.MELEME,MELEMA.MELEME
  37. POINTEUR IZLEMC.MELEME,MELEMS.MELEME,MELEMC.MELEME
  38. POINTEUR MELEMI.MELEME,MELEMP.MELEME
  39. POINTEUR MACRO1.MELEME
  40.  
  41. CHARACTER*8 TYPE,TYPC,NOMZ,NOMP,NOMD,NOM0
  42. CHARACTER*8 NOMPP,NOMDD
  43. CHARACTER*4 NOM
  44. INTEGER IPAD,IPAD2,IK
  45. REAL*8 KAUX,TETA1
  46. DIMENSION IXV(5)
  47. C
  48. C*************************************************************************
  49. CKMIC
  50. C WRITE(6,*)' Operateur KMIC MTABX=',MTABX
  51. C/MELEMS
  52.  
  53. C- Récupération de la table EQEX (pointeur MTAB1)
  54. C
  55. CALL LEKTAB(MTABX,'EQEX',MTAB1)
  56. IF(MTAB1.EQ.0)THEN
  57. C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
  58. MOTERR( 1: 8) = ' EQEX '
  59. MOTERR( 9:16) = ' EQEX '
  60. MOTERR(17:24) = ' KIZX '
  61. CALL ERREUR(786)
  62. RETURN
  63. ENDIF
  64. CALL ACME(MTAB1,'NAVISTOK',NASTOK)
  65. IF(NASTOK.EQ.0)THEN
  66. CALL ZKMIC(IKAS,MTABX,MTAB1)
  67. RETURN
  68. ENDIF
  69. C
  70. C- Récupération de la table INCO (pointeur KINC)
  71. C
  72. CALL LEKTAB(MTAB1,'INCO',KINC)
  73. IF(KINC.EQ.0)THEN
  74. C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
  75. MOTERR( 1: 8) = ' INCO '
  76. MOTERR( 9:16) = ' INCO '
  77. MOTERR(17:24) = ' EQEX '
  78. CALL ERREUR(786)
  79. RETURN
  80. ENDIF
  81.  
  82. C*************************************************************************
  83. C OPTIONS
  84. C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> CN
  85. C KFORM = 0 -> EFSI 1 -> EF 2 -> VF 3 -> EFMC
  86. C KPRE=3 pression P0 KPRE=4 pression P1 KPRE=2 cas macro 1ère génération
  87. C KPRE=5 pression MSOMMET
  88.  
  89. IAXI=0
  90. IK=0
  91. IF(IFOMOD.EQ.0)IAXI=2
  92. C
  93. C- Récupération de la table des options KOPT (pointeur KOPTI)
  94. C
  95. CALL LEKTAB(MTABX,'KOPT',KOPTI)
  96. IF (KOPTI.EQ.0) THEN
  97. C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
  98. MOTERR( 1: 8) = ' KOPT '
  99. MOTERR( 9:16) = ' KOPT '
  100. MOTERR(17:24) = ' KIZX '
  101. CALL ERREUR(786)
  102. RETURN
  103. ENDIF
  104.  
  105. CALL ACME(KOPTI,'KIMPL',KIMPL)
  106. CALL ACME(KOPTI,'KPOIN',KPRE)
  107. CALL ACMF(KOPTI,'AIMPL',TETA1)
  108. CALL ACME(KOPTI,'KFORM',KFORM)
  109. C gbm rajoute option décentrement 18/10/99
  110. CALL ACME(KOPTI,'IDCEN',IDCEN)
  111. C gbm rajoute IKOMP = 2 pour l'opérateur (div) et
  112. C on ne modifie pas IKOMP pour grad
  113. C pour appel ultérieur à SUPGEF.
  114. C IKOMP = 2 veut dire conservatif mais on a div(flux)
  115. C au lieu de div(ro*U).
  116. if (IKAS .EQ. 1 .OR. IKAS .EQ. 3) THEN
  117. IKOMP = 2
  118. endif
  119. * GBM modif 22/11/99
  120. CALL ACMF(KOPTI,'CMD',CMD)
  121.  
  122.  
  123. IF (IERR.NE.0) RETURN
  124. C WRITE(6,*)' Apres les options '
  125. C*************************************************************************
  126. C
  127. C- Récupération de la table DOMAINE associée au domaine local
  128. C
  129. CALL ACMM(MTABX,'NOMZONE',NOMZ)
  130. TYPE=' '
  131. CALL ACMO(MTABX,'DOMZ',TYPE,MMODEL)
  132. IF(TYPE.NE.'MMODEL')THEN
  133. C On attend un des objets : %m1:8 %m9:16 %m17:24 %m25:32 %m33:40
  134. MOTERR( 1: 8) = ' MMODEL '
  135. MOTERR( 8:16) = ' MMODEL '
  136. MOTERR(17:24) = ' MMODEL '
  137. MOTERR(25:32) = ' MMODEL '
  138. MOTERR(33:40) = ' MMODEL '
  139. CALL ERREUR(471)
  140. RETURN
  141. ENDIF
  142.  
  143. CALL LEKMOD(MMODEL,MTABZ,INEFMD)
  144.  
  145. CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
  146. CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
  147. CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
  148. C?? CALL LEKTAB(MTABZ,'MACRO',MACRO)
  149. MACRO1=0
  150. IF(INEFMD.EQ.2)CALL LEKTAB(MTABZ,'MACRO1',MACRO1)
  151. C?? IF(INEFMD.EQ.3)CALL LEKTAB(MTABZ,'QUADRATI',MQUAD)
  152. IF (IERR.NE.0) RETURN
  153. IF(INEFMD.EQ.4.AND.KPRE.NE.5)THEN
  154. C% Données incompatibles
  155. CALL ERREUR(21)
  156. RETURN
  157. ENDIF
  158.  
  159. MELEMI=MELEME
  160. IF(MACRO1.NE.0.AND.KPRE.NE.2)THEN
  161. MELEMI=MELEME
  162. ENDIF
  163.  
  164. IF(KPRE.EQ.2.AND.INEFMD.EQ.3)KPRE=3
  165. IF(INEFMD.EQ.1)KPRE=2
  166.  
  167. * ON MODIFIE POUR SOMMET -> SOMMET GBM 18/10/99
  168. MELEMP=MELEME
  169. CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
  170.  
  171. C*************************************************************************
  172. C VERIFICATIONS SUR LES INCONNUES
  173.  
  174. C WRITE(6,*)' Verification sur les inconnues '
  175. TYPE='LISTMOTS'
  176. CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO)
  177. IF (IERR.NE.0) RETURN
  178. SEGACT LINCO
  179. NBINC=LINCO.MOTS(/2)
  180. IF(NBINC.NE.2)THEN
  181. WRITE(IOIMP,*)'Operateur KMAC '
  182. WRITE(IOIMP,*)'Nombre d''inconnues incorrecte : ',NBINC,
  183. $ ' On attend 2'
  184. C Indice %m1:8 : contient plus de %i1 %m9:16
  185. MOTERR( 1:8) = 'LISTINCO'
  186. INTERR(1) = 2
  187. MOTERR(9:16) = ' MOTS '
  188. CALL ERREUR(799)
  189. RETURN
  190. ENDIF
  191.  
  192. C On recupere PHI n et TETA n pour Cranck-Nicholson
  193. NOMP=LINCO.MOTS(1)
  194. TYPE=' '
  195. CALL ACMO(KINC,NOMP,TYPE,MCHPOI)
  196. IF(TYPE.NE.'CHPOINT ')THEN
  197. WRITE(IOIMP,*)' Opérateur KMAC :'
  198. WRITE(IOIMP,*)' L objet CHPOINT ',NOMP,
  199. & ' n existe pas dans la table'
  200. C Indice %m1:8 : ne contient pas un objet de type %m9:16
  201. MOTERR( 1: 8) = 'INC '//NOMP
  202. MOTERR( 9:16) = 'CHPOINT '
  203. CALL ERREUR(800)
  204. RETURN
  205. ELSE
  206. CALL LICHT(MCHPOI,IZTU1,TYPC,IGEOM0)
  207. ENDIF
  208. C*************************************************************************
  209. C Le domaine de definition est donne par le SPG de la premiere inconnue
  210. C Les inconnues suivantes devront posseder ce meme pointeur
  211. C On verifie que les points de la zone sont tous inclus dans ce SPG
  212. C Inconnue Primale
  213.  
  214. C WRITE(6,*)' Verification inconnue primale '
  215. CALL KRIPAD(IGEOM0,IPADU)
  216. IF(IKAS.EQ.1.OR.IKAS.EQ.3)THEN
  217. MELEMK=MELEMS
  218. ELSE
  219. MELEMK=MELEMC
  220. ENDIF
  221.  
  222. CALL VERPAD(IPADU,MELEMK,IRET)
  223. IF(IRET.NE.0)THEN
  224. WRITE(IOIMP,*)' Opérateur KMAC '
  225. WRITE(IOIMP,*)' La zone ',NOMZ,
  226. $ ' n''est pas incluse dans le domaine'
  227. & , ' de définition de l''inconnue ',NOMP
  228. WRITE(IOIMP,*)' MELEMK=',melemk,' IGEOM0=',IGEOM0
  229. C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
  230. MOTERR(1: 8) = 'INC '//NOMP
  231. MOTERR(9:16) = 'CHPOINT '
  232. CALL ERREUR(788)
  233. IPAS=0
  234. RETURN
  235. ENDIF
  236.  
  237.  
  238. C*************************************************************************
  239.  
  240. NOMD=LINCO.MOTS(2)
  241. TYPE=' '
  242. CALL ACMO(KINC,NOMD,TYPE,MCHPOI)
  243. IF(TYPE.NE.'CHPOINT ')THEN
  244. WRITE(IOIMP,*)' Opérateur KMAC :'
  245. WRITE(IOIMP,*)' L objet CHPOINT ',NOMD,
  246. & ' n existe pas dans la table'
  247. C Indice %m1:8 : ne contient pas un objet de type %m9:16
  248. MOTERR( 1: 8) = 'INC '//NOMD
  249. MOTERR( 9:16) = 'CHPOINT '
  250. CALL ERREUR(800)
  251. RETURN
  252. ELSE
  253. CALL LICHT(MCHPOI,TETAN,TYPC,IGEOM0)
  254. ENDIF
  255.  
  256. NC=TETAN.VPOCHA(/2)
  257. C*************************************************************************
  258. C Le domaine de definition est donne par le SPG de la premiere inconnue
  259. C Les inconnues suivantes devront posseder ce meme pointeur
  260. C On verifie que les points de la zone sont tous inclus dans ce SPG
  261. C Inconnue Duale
  262.  
  263. C WRITE(6,*)' IGEOM0=',igeom0
  264. CALL KRIPAD(IGEOM0,MLENTI)
  265. IF(IKAS.EQ.1.OR.IKAS.EQ.3)THEN
  266. C ON MODIFIE POUR CAS MEME ESPACE DIV U SOMMET -> SOMMET GBM 18/10
  267. MELEMK=MELEMS
  268. ELSE
  269. MELEMK=MELEMK
  270. ENDIF
  271.  
  272. C WRITE(6,*)' Verification inconnue duale ',MELEMK
  273. CALL VERPAD(MLENTI,MELEMK,IRET)
  274. IF(IRET.NE.0)THEN
  275. WRITE(IOIMP,*)' Opérateur KMAC '
  276. WRITE(IOIMP,*)' La zone ',NOMZ,
  277. $ ' n''est pas incluse dans le domaine'
  278. & ,' de définition de l''inconnue ',NOMD
  279. WRITE(IOIMP,*)' MELEMK=',melemk,' IGEOM0=',IGEOM0
  280. C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
  281. MOTERR(1: 8) = 'INC '//NOMD
  282. MOTERR(9:16) = 'CHPOINT '
  283. CALL ERREUR(788)
  284. IPAS=0
  285. RETURN
  286. ENDIF
  287.  
  288. C GBM rajoute NPTU 18/10
  289. NPTU = IZTU1.VPOCHA(/1)
  290.  
  291. C*************************************************************************
  292. C Lecture du ou des coefficients
  293. C Type du coefficient :
  294. C IK1=0 CHPOINT IK1=1 scalaire IK1=2 vecteur
  295.  
  296. C WRITE(6,*)' Verification sur les coefficients '
  297. CALL ACME(MTABX,'IARG',IARG)
  298. C WRITE(6,*)' IARG=',iarg
  299.  
  300. IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)THEN
  301. TYPE=' '
  302. CALL ACMO(MTABZ,'MELSTB',TYPE,MELSTB)
  303. SEGACT MELSTB
  304. IF(IDIM.EQ.2)NBELEM=MELSTB.NUM(/2)/4
  305. IF(IDIM.EQ.3)NBELEM=MELSTB.NUM(/2)/8
  306. NBNN=MELSTB.NUM(/1)
  307. NBSOUS=0
  308. NBREF=0
  309. SEGINI MELEMA
  310. MELEMA.ITYPEL=MELSTB.ITYPEL
  311.  
  312. NKPE=4
  313. IF(IDIM.EQ.3)NKPE=8
  314. do 4878 k=1,nbelem
  315. mi=(k-1)*NKPE+1
  316. do 4879 i=1,nbnn
  317. MELEMA.num(i,k)=melstb.num(i,mi)
  318. 4879 continue
  319. C WRITE(6,*)k,(MELEMA.num(i,k),i=1,nbnn)
  320. 4878 continue
  321.  
  322. TYPE=' '
  323. CALL ACMO(MTABZ,'MCHPOC',TYPE,MCHPOC)
  324. TYPE=' '
  325. CALL ACMO(MTABZ,'CENTRE',TYPE,MCTREI)
  326. ENDIF
  327.  
  328.  
  329. C 1er COEF
  330.  
  331. IXV(1)=MELEMC
  332. IXV(2)=1
  333. IXV(3)=0
  334. IXV(4)=MELEMS
  335. IXV(5)=-MELEMS
  336. IRET=5
  337. CALL LEKCOF('Opérateur KMAC :',
  338. & MTABX,KINC,1,IXV,IZTG1,IZTGG1,NPT1,NC1,IK1,IRET)
  339. IF(IRET.EQ.0)RETURN
  340. IF(IK1.GE.4)CALL KRIPAD(MELEMS,MLENTI)
  341.  
  342.  
  343. BETA0=1.D-1
  344. IF(IARG.EQ.2)THEN
  345. C 2ème COEF
  346. IXV(1)=0
  347. IXV(2)=1
  348. IXV(3)=0
  349. CALL LEKCOF('Opérateur KMAC :',
  350. & MTABX,KINC,2,IXV,IZTG2,IZBETA,NPT2,NC2,IK2,IRET)
  351. IF(IRET.EQ.0)RETURN
  352. BETA0=IZBETA.VPOCHA(1,1)
  353. ENDIF
  354.  
  355. C GBM rajoute 3eme coef mot clef 'SOMM'
  356.  
  357. C GBM rajoute 4eme coef viscosité pour SUPGEF
  358.  
  359. ANUK=1.D-15
  360. IF(IARG.EQ.4)THEN
  361. C 4ème COEF
  362. IXV(1)=0
  363. IXV(2)=1
  364. IXV(3)=0
  365. CALL LEKCOF('Opérateur KMAC :',
  366. & MTABX,KINC,4,IXV,IZTGG,IZANU,NPT3,NC3,IK3,IRET)
  367. IF(IRET.EQ.0)RETURN
  368. ANUK=IZANU.VPOCHA(1,1)
  369. ENDIF
  370.  
  371. NOMP=LINCO.MOTS(1)
  372. NOMD=LINCO.MOTS(2)(1:4)
  373.  
  374. NRIGE=7
  375. NKID =9
  376. NKMT =7
  377. NMATRI=1
  378. IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)NMATRI=2
  379. IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.4.AND.IDIM.EQ.2)NMATRI=2
  380. SEGINI MATRIK
  381.  
  382. C Fin Stabilisation de toutes sortes
  383. C WRITE(6,*)'C Fin Stabilisation de toutes sortes'
  384.  
  385. NBME=IDIM
  386. CALL KRIPAD(MELEMI,MLENTI1)
  387. SEGACT MELEMI
  388. NBSOUS=MELEMI.LISOUS(/1)
  389. IF(NBSOUS.EQ.0)NBSOUS=1
  390. SEGINI IMATRI
  391.  
  392. IF(IKAS.EQ.2)THEN
  393. KSPGD=MELEMS
  394. KSPGP=MELEMC
  395. IRIGEL(2,1)=MELEME
  396. IRIGEL(1,1)=MELEME
  397. ELSE
  398. KSPGP=MELEMS
  399. KSPGD=MELEMS
  400. IRIGEL(1,1)=MELEME
  401. IRIGEL(2,1)=MELEME
  402. ENDIF
  403. SEGACT MELEMP
  404.  
  405. C WRITE(6,*)' ds kmac melemp=',IRIGEL(1,1)
  406. C WRITE(6,*)' ds kmac melemd=',IRIGEL(2,1)
  407.  
  408. IRIGEL(4,1)=IMATRI
  409. * je modifie celui là pour retirer les multiplicateurs
  410. IF(IKAS.EQ.1)IRIGEL(7,1)= 5
  411. IF(IKAS.EQ.2)IRIGEL(7,1)=5
  412. IF(IKAS.EQ.3)IRIGEL(7,1)=5
  413.  
  414. K0=0
  415. DO 11 L=1,NBSOUS
  416. IPT1=MELEMI
  417. IF(NBSOUS.NE.1)IPT1=MELEMI.LISOUS(L)
  418. SEGACT IPT1
  419. NBEL=IPT1.NUM(/2)
  420.  
  421. IF(IKAS.EQ.2)THEN
  422. MP=IPT1.NUM(/1)
  423. NP=MELEMP.NUM(/1)
  424. ELSE
  425. NP=IPT1.NUM(/1)
  426. MP=MELEMP.NUM(/1)
  427. ENDIF
  428.  
  429. DO 12 I=1,NBME
  430. SEGINI IZAFM
  431. C WRITE(6,*)' ni izafm np=',np,' mp=',mp,' nbel=',nbel,izafm,l,i
  432. LIZAFM(L,I)=IZAFM
  433. IF(IKAS.EQ.2)THEN
  434. WRITE(NOM,FMT='(I1,A3)')I,NOMD(1:3)
  435. LISDUA(I)=NOM//' '
  436. LISPRI(I)=NOMP
  437. ELSE
  438. WRITE(NOM,FMT='(I1,A3)')I,NOMP(1:3)
  439. LISPRI(I)=NOM//' '
  440. LISDUA(I)=NOMD
  441. ENDIF
  442. 12 CONTINUE
  443. IPM1=LIZAFM(L,1)
  444. IPM2=LIZAFM(L,2)
  445. IPM3=LIZAFM(L,2)
  446. IF(IDIM.EQ.3)IPM3=LIZAFM(L,3)
  447.  
  448. IF (IKAS .EQ. 1 .OR. IKAS .EQ. 3) THEN
  449. c on code div roU -> T,
  450. IF(IK1.LT.4)THEN
  451. c cas où le coef multiplicateur est scalaire, point
  452. c ou champ au centre.
  453. CALL KSPRUS
  454. & (IPT1,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,IK1,K0,
  455. & ANUK,IDCEN,IKOMP,NPTU,MLENTI,IZTU1,TETAN,NP,NBEL,CMD)
  456. ELSE
  457. C MODIFé par GBM, cas coef au sommet.
  458. CALL KSPRJS
  459. & (IPT1,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1
  460. & ,IPADU,IK1)
  461. ENDIF
  462. ELSE
  463. c on code gradT -> U
  464. c on inverse la place de IZTU1 et de TETAN dans les arguments
  465. c ca l'utilisateur rentre la pression en premier, IZTU1 est donc
  466. c une pression (nom mal choisi). GBM 18/10/99
  467. IF(IK1.LT.4)THEN
  468. c cas où le coef multiplicateur est scalaire, point
  469. c ou champ au centre.
  470. CALL KSPRUS
  471. & (IPT1,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,IK1,K0,
  472. & ANUK,IDCEN,IKOMP,NPTU,MLENTI,TETAN,IZTU1,NP,NBEL,CMD)
  473. ELSE
  474. C MODIFé par GBM, cas coef au sommet.
  475. CALL KSPRJS
  476. & (IPT1,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,
  477. & IPADU,IK1)
  478. ENDIF
  479. ENDIF
  480. C
  481. K0=K0+NBEL
  482. SEGDES IPM1,IPM2,IPM3
  483. SEGDES IPT1
  484. 11 CONTINUE
  485. SEGDES MELEMI
  486. IF(IK1.EQ.4)SEGSUP MLENTI
  487.  
  488. CALL ECROBJ('MATRIK',MATRIK)
  489.  
  490. NAT=2
  491. NSOUPO=0
  492. SEGINI MCHPOI
  493. JATTRI(1)=2
  494. CALL ECROBJ('CHPOINT',MCHPOI)
  495.  
  496. C WRITE(6,*)' Fin operateur KMIC'
  497. SEGDES IMATRI,MATRIK
  498. RETURN
  499. 1001 FORMAT(20(1X,I5))
  500. END
  501.  
  502.  
  503.  
  504.  
  505.  
  506.  
  507.  
  508.  
  509.  
  510.  
  511.  
  512.  
  513.  
  514.  
  515.  
  516.  
  517.  
  518.  

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