Télécharger divs.eso

Retour à la liste

Numérotation des lignes :

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

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