Télécharger kmic.eso

Retour à la liste

Numérotation des lignes :

  1. C KMIC SOURCE PV 16/11/17 21:59:48 9180
  2. SUBROUTINE KMIC(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,VISCO.MPOVAL
  27. POINTEUR IPM.IZAFM
  28. SEGMENT IMATRS
  29. ENDSEGMENT
  30. POINTEUR IPMS.IZAFM,IPS1.IZAFM,IPS2.IZAFM,IPS3.IZAFM
  31.  
  32. PARAMETER (LRV=64)
  33. DIMENSION VISC(3)
  34.  
  35. -INC SMLENTI
  36. POINTEUR IZIPAD.MLENTI,MLENTI1.MLENTI,MLENTI2.MLENTI
  37. -INC SMLMOTS
  38. POINTEUR LINCO.MLMOTS
  39. -INC SMELEME
  40. POINTEUR MELEMZ.MELEME,MELEMB.MELEME,MELSTB.MELEME
  41. POINTEUR MELEM1.MELEME,MELES1.MELEME,MCTREI.MELEME
  42. POINTEUR IGEOM.MELEME,MELEMA.MELEME
  43. POINTEUR IZLEMC.MELEME,MELEMS.MELEME,MELEMC.MELEME
  44. POINTEUR MELEMI.MELEME,MELEMP.MELEME
  45. POINTEUR MACRO1.MELEME
  46.  
  47. CHARACTER*8 TYPE,TYPC,NOMZ,NOMP,NOMD,NOM0,NOMK
  48. CHARACTER*8 NOMPP,NOMDD
  49. CHARACTER*4 NOM
  50. INTEGER IPAD,IPAD2,IK
  51. REAL*8 KAUX,TETA1
  52. DIMENSION IXV(5)
  53. C
  54. C*************************************************************************
  55. CKMIC
  56. C write(6,*)' Operateur KMIC MTABX=',MTABX
  57. C
  58. C- Récupération de la table EQEX (pointeur MTAB1)
  59. C
  60. CALL LEKTAB(MTABX,'EQEX',MTAB1)
  61. IF(MTAB1.EQ.0)THEN
  62. C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
  63. MOTERR( 1: 8) = ' EQEX '
  64. MOTERR( 9:16) = ' EQEX '
  65. MOTERR(17:24) = ' KIZX '
  66. CALL ERREUR(786)
  67. RETURN
  68. ENDIF
  69.  
  70. CALL ACME(MTABX,'IARG',IARG)
  71. IF(IARG.GE.3)THEN
  72. C 3ème COEF
  73. CALL ACMM(MTABX,'ARG3',NOMK)
  74.  
  75. IF(NOMK.EQ.'SOMM')THEN
  76. * Cet opérateur calcul un gradient ou div de sommets vers sommets
  77. * sans contrainte de type multiplicateur de Lagrange.
  78. CALL DIVS(IKAS,MTABX)
  79. RETURN
  80. *
  81. ELSEIF(NOMK.EQ.'NOMU')THEN
  82. * Cet opérateur calcule la divergence (ou grad) de sommets vers ctp1
  83. * sans contrainte de cmultiplicateurs de Lagrange.
  84. CALL DIVC(IKAS,MTABX)
  85. C write(6,*)' RETOUR KKMIC'
  86. RETURN
  87.  
  88. ENDIF
  89. ENDIF
  90.  
  91. CALL ACME(MTAB1,'DISCPRES',IDISCP)
  92.  
  93. CALL ACME(MTAB1,'NAVISTOK',NASTOK)
  94. IF(NASTOK.EQ.0)THEN
  95. CALL ZKMIC(IKAS,MTABX,MTAB1)
  96. RETURN
  97. ENDIF
  98. C
  99. C- Récupération de la table INCO (pointeur KINC)
  100. C
  101. CALL LEKTAB(MTAB1,'INCO',KINC)
  102. IF(KINC.EQ.0)THEN
  103. C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
  104. MOTERR( 1: 8) = ' INCO '
  105. MOTERR( 9:16) = ' INCO '
  106. MOTERR(17:24) = ' EQEX '
  107. CALL ERREUR(786)
  108. RETURN
  109. ENDIF
  110.  
  111. C*************************************************************************
  112. C OPTIONS
  113. C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> CN
  114. C KFORM = 0 -> EFSI 1 -> EF 2 -> VF 3 -> EFMC
  115. C KPRE=3 pression P0 KPRE=4 pression P1 KPRE=2 cas macro 1ère génération
  116. C KPRE=5 pression MSOMMET
  117.  
  118. NOMUL=1
  119. IAXI=0
  120. IK=0
  121. IF(IFOMOD.EQ.0)IAXI=2
  122. C
  123. C- Récupération de la table des options KOPT (pointeur KOPTI)
  124. C
  125. CALL LEKTAB(MTABX,'KOPT',KOPTI)
  126. IF (KOPTI.EQ.0) THEN
  127. C Indice %m1:8 : Indice %m9:16 non trouvé dans la table %m17:24
  128. MOTERR( 1: 8) = ' KOPT '
  129. MOTERR( 9:16) = ' KOPT '
  130. MOTERR(17:24) = ' KIZX '
  131. CALL ERREUR(786)
  132. RETURN
  133. ENDIF
  134.  
  135. CALL ACME(KOPTI,'KIMPL',KIMPL)
  136. CALL ACME(KOPTI,'KPOIN',KPRE)
  137. C write(6,*)'KBBT KPRE=',KPRE
  138. CALL ACMF(KOPTI,'AIMPL',TETA1)
  139. CALL ACME(KOPTI,'KFORM',KFORM)
  140. C CALL ACME(KOPTI,'IKOMP',IKOMP)
  141.  
  142. IF (IERR.NE.0) RETURN
  143. C write(6,*)' Apres les options '
  144. C*************************************************************************
  145. C
  146. C- Récupération de la table DOMAINE associée au domaine local
  147. C
  148. CALL ACMM(MTABX,'NOMZONE',NOMZ)
  149. TYPE=' '
  150. CALL ACMO(MTABX,'DOMZ',TYPE,MMODEL)
  151. IF(TYPE.NE.'MMODEL')THEN
  152. C On attend un des objets : %m1:8 %m9:16 %m17:24 %m25:32 %m33:40
  153. MOTERR( 1: 8) = ' MMODEL '
  154. MOTERR( 8:16) = ' MMODEL '
  155. MOTERR(17:24) = ' MMODEL '
  156. MOTERR(25:32) = ' MMODEL '
  157. MOTERR(33:40) = ' MMODEL '
  158. CALL ERREUR(471)
  159. RETURN
  160. ENDIF
  161.  
  162. CALL LEKMOD(MMODEL,MTABZ,INEFMD)
  163.  
  164. CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
  165. CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
  166. CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
  167. C?? CALL LEKTAB(MTABZ,'MACRO',MACRO)
  168. MACRO1=0
  169. IF(INEFMD.EQ.2)CALL LEKTAB(MTABZ,'MACRO1',MACRO1)
  170. C?? IF(INEFMD.EQ.3)CALL LEKTAB(MTABZ,'QUADRATI',MQUAD)
  171. IF (IERR.NE.0) RETURN
  172. IF(INEFMD.EQ.4.AND.KPRE.NE.5)THEN
  173. C% Données incompatibles
  174. CALL ERREUR(21)
  175. RETURN
  176. ENDIF
  177.  
  178. MELEMI=MELEME
  179. IF(MACRO1.NE.0.AND.KPRE.NE.2)THEN
  180. MELEMI=MACRO1
  181. ENDIF
  182.  
  183. IF(KPRE.EQ.2.AND.INEFMD.EQ.3)KPRE=3
  184. IF(INEFMD.EQ.1.AND.KPRE.NE.5)KPRE=2
  185.  
  186. IF(KPRE.EQ.3)THEN
  187. CALL LEKTAB(MTABZ,'CENTREP0',MELEMC)
  188. MELEMP=MELEMC
  189. ELSEIF(KPRE.EQ.4)THEN
  190. CALL LEKTAB(MTABZ,'CENTREP1',MELEMC)
  191. CALL LEKTAB(MTABZ,'ELTP1NC ',MELEMP)
  192. ELSEIF(KPRE.EQ.2)THEN
  193. CALL LEKTAB(MTABZ,'CENTRE',MELEMC)
  194. MELEMP=MELEMC
  195. ELSEIF(KPRE.EQ.5)THEN
  196. CALL LEKTAB(MTABZ,'MSOMMET',MELEMC)
  197. CALL LEKTAB(MTABZ,'MMAIL ',MELEMP)
  198. ENDIF
  199.  
  200. IF(IDISCP.EQ.0)THEN
  201. CALL ECME(MTAB1,'DISCPRES',KPRE)
  202. ELSEIF(IDISCP.NE.KPRE)THEN
  203. C% Données incompatibles
  204. CALL ERREUR(21)
  205. RETURN
  206. ENDIF
  207.  
  208. IKOMP=1
  209. C IF(KPRE.EQ.5)IKOMP=0
  210. C IF(KPRE.EQ.5)IKOMP=1
  211.  
  212. C*************************************************************************
  213. C VERIFICATIONS SUR LES INCONNUES
  214.  
  215. C write(6,*)' Verification sur les inconnues '
  216. TYPE='LISTMOTS'
  217. CALL ACMO(MTABX,'LISTINCO',TYPE,LINCO)
  218. IF (IERR.NE.0) RETURN
  219. SEGACT LINCO
  220. NBINC=LINCO.MOTS(/2)
  221. IF(NBINC.NE.2)THEN
  222. WRITE(6,*)'Operateur KMAB '
  223. WRITE(6,*)'Nombre d''inconnues incorrecte : ',NBINC,' On attend 2'
  224. C Indice %m1:8 : contient plus de %i1 %m9:16
  225. MOTERR( 1:8) = 'LISTINCO'
  226. INTERR(1) = 2
  227. MOTERR(9:16) = ' MOTS '
  228. CALL ERREUR(799)
  229. RETURN
  230. ENDIF
  231.  
  232. C On recupere PHI n et TETA n pour Cranck-Nicholson
  233. NOMP=LINCO.MOTS(1)
  234. TYPE=' '
  235. CALL ACMO(KINC,NOMP,TYPE,MCHPOI)
  236. IF(TYPE.NE.'CHPOINT ')THEN
  237. WRITE(6,*)' Opérateur KMAB :'
  238. WRITE(6,*)' L objet CHPOINT ',NOMP,
  239. & ' n existe pas dans la table'
  240. C Indice %m1:8 : ne contient pas un objet de type %m9:16
  241. MOTERR( 1: 8) = 'INC '//NOMP
  242. MOTERR( 9:16) = 'CHPOINT '
  243. CALL ERREUR(800)
  244. RETURN
  245. ELSE
  246. CALL LICHT(MCHPOI,IZTU1,TYPC,IGEOM0)
  247. ENDIF
  248. C*************************************************************************
  249. C Le domaine de definition est donne par le SPG de la premiere inconnue
  250. C Les inconnues suivantes devront posseder ce meme pointeur
  251. C On verifie que les points de la zone sont tous inclus dans ce SPG
  252. C Inconnue Primale
  253.  
  254. C write(6,*)' Verification inconnue primale '
  255. CALL KRIPAD(IGEOM0,MLENTI)
  256. IF(IKAS.EQ.1.OR.IKAS.EQ.3)THEN
  257. MELEMK=MELEMS
  258. ELSE
  259. MELEMK=MELEMC
  260. ENDIF
  261.  
  262. C write(6,*)' MELEMK=',melemk
  263. CALL VERPAD(MLENTI,MELEMK,IRET)
  264. IF(IRET.NE.0)THEN
  265. WRITE(6,*)' Opérateur KMAB '
  266. WRITE(6,*)' La zone ',NOMZ,' n''est pas incluse dans le domaine'
  267. & , ' de définition de l''inconnue ',NOMP
  268. WRITE(6,*)' MELEMK=',melemk,' IGEOM0=',IGEOM0
  269. C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
  270. MOTERR(1: 8) = 'INC '//NOMP
  271. MOTERR(9:16) = 'CHPOINT '
  272. CALL ERREUR(788)
  273. IPAS=0
  274. RETURN
  275. ENDIF
  276.  
  277. SEGSUP MLENTI
  278.  
  279. C*************************************************************************
  280.  
  281. NOMD=LINCO.MOTS(2)
  282. TYPE=' '
  283. CALL ACMO(KINC,NOMD,TYPE,MCHPOI)
  284. IF(TYPE.NE.'CHPOINT ')THEN
  285. WRITE(6,*)' Opérateur KMAB :'
  286. WRITE(6,*)' L objet CHPOINT ',NOMD,
  287. & ' n existe pas dans la table'
  288. C Indice %m1:8 : ne contient pas un objet de type %m9:16
  289. MOTERR( 1: 8) = 'INC '//NOMD
  290. MOTERR( 9:16) = 'CHPOINT '
  291. CALL ERREUR(800)
  292. RETURN
  293. ELSE
  294. CALL LICHT(MCHPOI,TETAN,TYPC,IGEOM0)
  295. ENDIF
  296.  
  297. NC=TETAN.VPOCHA(/2)
  298. C*************************************************************************
  299. C Le domaine de definition est donne par le SPG de la premiere inconnue
  300. C Les inconnues suivantes devront posseder ce meme pointeur
  301. C On verifie que les points de la zone sont tous inclus dans ce SPG
  302. C Inconnue Duale
  303.  
  304. C write(6,*)' IGEOM0=',igeom0
  305. CALL KRIPAD(IGEOM0,MLENTI)
  306. IF(IKAS.EQ.1.OR.IKAS.EQ.3)THEN
  307. MELEMK=MELEMC
  308. ELSE
  309. MELEMK=MELEMS
  310. ENDIF
  311.  
  312. C write(6,*)' Verification inconnue duale ',MELEMK
  313. CALL VERPAD(MLENTI,MELEMK,IRET)
  314. IF(IRET.NE.0)THEN
  315. WRITE(6,*)' Opérateur KMAB '
  316. WRITE(6,*)' La zone ',NOMZ,' n''est pas incluse dans le domaine'
  317. & ,' de définition de l''inconnue ',NOMD
  318. WRITE(6,*)' MELEMK=',melemk,' IGEOM0=',IGEOM0
  319. C Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
  320. MOTERR(1: 8) = 'INC '//NOMD
  321. MOTERR(9:16) = 'CHPOINT '
  322. CALL ERREUR(788)
  323. IPAS=0
  324. RETURN
  325. ENDIF
  326.  
  327. SEGSUP MLENTI
  328.  
  329. C*************************************************************************
  330. C Lecture du ou des coefficients
  331. C Type du coefficient :
  332. C IK1=0 CHPOINT IK1=1 scalaire IK1=2 vecteur
  333.  
  334. C write(6,*)' Verification sur les coefficients '
  335. CALL ACME(MTABX,'IARG',IARG)
  336. C write(6,*)' IARG=',iarg
  337.  
  338. IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)THEN
  339. TYPE=' '
  340. CALL ACMO(MTABZ,'MELSTB',TYPE,MELSTB)
  341. SEGACT MELSTB
  342. IF(IDIM.EQ.2)NBELEM=MELSTB.NUM(/2)/4
  343. IF(IDIM.EQ.3)NBELEM=MELSTB.NUM(/2)/8
  344. NBNN=MELSTB.NUM(/1)
  345. NBSOUS=0
  346. NBREF=0
  347. SEGINI MELEMA
  348. MELEMA.ITYPEL=MELSTB.ITYPEL
  349.  
  350. NKPE=4
  351. IF(IDIM.EQ.3)NKPE=8
  352. do 4878 k=1,nbelem
  353. mi=(k-1)*NKPE+1
  354. do 4879 i=1,nbnn
  355. MELEMA.num(i,k)=melstb.num(i,mi)
  356. 4879 continue
  357. C write(6,*)k,(MELEMA.num(i,k),i=1,nbnn)
  358. 4878 continue
  359.  
  360. TYPE=' '
  361. CALL ACMO(MTABZ,'MCHPOC',TYPE,MCHPOC)
  362. TYPE=' '
  363. CALL ACMO(MTABZ,'CENTRE',TYPE,MCTREI)
  364. ENDIF
  365.  
  366.  
  367. C 1er COEF
  368.  
  369. IXV(1)=MELEMC
  370. IXV(2)=1
  371. IXV(3)=0
  372. IXV(4)=MELEMS
  373. IXV(5)=-MELEMS
  374. IRET=5
  375. CALL LEKCOF('Opérateur KMAB :',
  376. & MTABX,KINC,1,IXV,IZTG1,IZTGG1,NPT1,NC1,IK1,IRET)
  377. IF(IRET.EQ.0)RETURN
  378. IF(IK1.GE.4)CALL KRIPAD(MELEMS,MLENTI)
  379.  
  380.  
  381. BETA0=1.D-1
  382. IF(IARG.EQ.2)THEN
  383. C 2ème COEF
  384. IXV(1)=0
  385. IXV(2)=1
  386. IXV(3)=0
  387. CALL LEKCOF('Opérateur KMAB :',
  388. & MTABX,KINC,2,IXV,IZTG2,IZBETA,NPT2,NC2,IK2,IRET)
  389. IF(IRET.EQ.0)RETURN
  390. BETA0=IZBETA.VPOCHA(1,1)
  391. C? IF(IZBETA.VPOCHA(1,1).LT.0.D0)THEN
  392. C% Données incompatibles
  393. C? CALL ERREUR(21)
  394. C? RETURN
  395. C? ENDIF
  396. ENDIF
  397.  
  398. IF(IARG.EQ.3)THEN
  399. C 3ème COEF
  400. CALL ACMM(MTABX,'ARG3',NOMK)
  401. C write(6,*)' NOMK=',nomk
  402. IF(NOMK.EQ.'CAVIFERM')IKOMP=0
  403. IF(NOMK.EQ.'NOMUL ')NOMUL=0
  404. ENDIF
  405.  
  406. NOMP=LINCO.MOTS(1)
  407. NOMD=LINCO.MOTS(2)(1:4)
  408.  
  409. NRIGE=7
  410. NKID =9
  411. NKMT =7
  412. NMATRI=1
  413. IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)NMATRI=2
  414. IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.4.AND.IDIM.EQ.2)NMATRI=2
  415. IF(KPRE.EQ.5)NMATRI=NMATRI+1
  416. C write(6,*)' NMATRI=',nmatri,' KPRE=',kpre
  417. SEGINI MATRIK
  418.  
  419. C CAS Stabilisation via MACRO CENTRE
  420. IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.2)THEN
  421. IF(IKAS.EQ.3)BETA0=-BETA0
  422. C write(6,*)' CAS Stabilisation via MACRO CENTRE ',BETA0,IKAS
  423. NK=0
  424. I2=2
  425. NBME=1
  426. NBSOUS=1
  427. SEGINI IMATRI
  428. IRIGEL(4,i2)=IMATRI
  429. KSPGP=MCTREI
  430. KSPGD=MCTREI
  431. IRIGEL(1,i2)=MELEMA
  432. IRIGEL(2,i2)=MELEMA
  433. IRIGEL(7,i2)=0
  434. CALL LICHT(MCHPOC,MPOVAL,TYPC,IGEOM)
  435.  
  436. C write(6,*)' MELSTB=',melstb,' MCTREI=',MCTREI
  437.  
  438. SEGACT MELSTB
  439.  
  440. NBSOUS=MELSTB.LISOUS(/1)
  441. IF(NBSOUS.NE.0)THEN
  442. CALL ERREUR(5)
  443. ENDIF
  444.  
  445. NBEL=MELEMA.NUM(/2)
  446. NBCI=MELSTB.NUM(/2)
  447. NP =MELSTB.NUM(/1)
  448. MP =NP
  449.  
  450. C write(6,*)' nbel,nbci,np,mp=',nbel,nbci,np,mp
  451. SEGINI IZAFM
  452. LIZAFM(1,1)=IZAFM
  453. LISPRI(1)=NOMD
  454. LISDUA(1)=NOMD
  455.  
  456. CALL KRIPAD(MCTREI,MLENT1)
  457.  
  458. C write(6,*)' IK1=',ik1
  459. IF(IK1.LT.4)THEN
  460. DO 33 K=1,NBEL
  461.  
  462. NK=NK+1
  463. KC=1+(1-IK1)*(NK-1)
  464. C write(6,*)' K=',K,' NK=',nk,np,kc
  465.  
  466. DO 32 J=1,NP
  467. K1=MLENT1.LECT(MELEMA.NUM(J,K))
  468. II=J
  469. DO 321 I=1,NP
  470. U=VPOCHA(K1,I)*BETA0*IZTGG1.VPOCHA(KC,1)
  471. IF(I.EQ.1)U=ABS(VPOCHA(K1,I))*BETA0*IZTGG1.VPOCHA(KC,1)
  472. IF(II.LE.NP)THEN
  473. AM(K,II,J)=U
  474. ELSE
  475. AM(K,II-NP,J)=U
  476. ENDIF
  477. II=II+1
  478. 321 CONTINUE
  479. 32 CONTINUE
  480. 33 CONTINUE
  481. C write(6,*)' FIN BCL 33 '
  482. ELSE
  483.  
  484. DO 34 K=1,NBEL
  485.  
  486. NK=NK+1
  487. UA=1.D0
  488.  
  489. DO 35 J=1,NP
  490. K1=MLENT1.LECT(MELEMA.NUM(J,K))
  491. II=J
  492. DO 351 I=1,NP
  493. U=VPOCHA(K1,I)*BETA0*UA
  494. IF(I.EQ.1)U=ABS(VPOCHA(K1,I))*BETA0*UA
  495. IF(II.LE.NP)THEN
  496. AM(K,II,J)=U
  497. ELSE
  498. AM(K,II-NP,J)=U
  499. ENDIF
  500. II=II+1
  501. 351 CONTINUE
  502. 35 CONTINUE
  503. 34 CONTINUE
  504.  
  505.  
  506. ENDIF
  507. SEGSUP MLENT1
  508. ENDIF
  509. C write(6,*)' Apres 34 '
  510.  
  511. C CAS Stabilisation via MACRO CENTREP1
  512. DEUPI=1.D0
  513. IF(IAXI.NE.0)DEUPI=2.D0*XPI
  514. IF(IKAS.EQ.3)DEUPI=-DEUPI
  515.  
  516. IF(MACRO1.NE.0.AND.IKAS.NE.2.AND.KPRE.EQ.4.AND.IDIM.EQ.2)THEN
  517. I2=2
  518. NBME=1
  519. NBSOUS=1
  520. SEGINI IMATRI
  521. IRIGEL(4,i2)=IMATRI
  522. KSPGP=MELEMC
  523. KSPGD=MELEMC
  524. IRIGEL(1,i2)=MELEMP
  525. IRIGEL(2,i2)=MELEMP
  526. IRIGEL(7,i2)=0
  527.  
  528. SEGACT MELEMP
  529. NBSOUS=MELEMP.LISOUS(/1)
  530. C write(6,*)' MELEMP2 NBSOUS=',NBSOUS,' MELEMP=',MELEMP
  531. IF(NBSOUS.NE.0)THEN
  532. CALL ERREUR(5)
  533. ENDIF
  534.  
  535. SEGACT MELEMI
  536.  
  537. NBEL=MELEMP.NUM(/2)
  538. NP =MELEMP.NUM(/1)
  539. MP =NP
  540. NOM0=NOMS(MELEMP.ITYPEL)//' '
  541.  
  542. SEGINI IZAFM
  543. LIZAFM(1,1)=IZAFM
  544. LISPRI(1)=NOMD
  545. LISDUA(1)=NOMD
  546.  
  547. CALL KALPBG(NOM0,'FONFORM ',IZFFM)
  548.  
  549. SEGACT IZFFM*MOD
  550. IZHR=KZHR(1)
  551. SEGACT IZHR*MOD
  552. NES=GR(/1)
  553. NPG=GR(/3)
  554.  
  555. NK=0
  556. DO 430 LM=1,MAX(1,MELEMI.LISOUS(/1))
  557. IPT1=MELEMI
  558. IF(MELEMI.LISOUS(/1).NE.0)IPT1=MELEMI.LISOUS(LM)
  559. SEGACT IPT1
  560. ITYP=IPT1.ITYPEL
  561. BETA=0.D0
  562. IF(NOMS(ITYP).EQ.'TRI6')BETA=BETA0
  563. C write(6,*)' NOMS(ITYP)=',NOMS(ITYP),BETA
  564. NBEL=IPT1.NUM(/2)
  565.  
  566. C Cas coefficient FLOTTANT ou CENTRE
  567. IF(IK1.LT.4)THEN
  568. DO 43 KK=1,NBEL
  569. NK=NK+1
  570. DO 42 I=1,NP
  571. J=MELEMP.NUM(I,NK)
  572. DO 42 N=1,IDIM
  573. XYZ(N,I)=XCOOR((J-1)*(IDIM+1) +N)
  574. 42 CONTINUE
  575.  
  576. CALL CALJBR
  577. &(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
  578. KC=1+(1-IK1)*(NK-1)
  579.  
  580. DO 44 I=1,NP
  581. DO 44 J=1,NP
  582. U=0.D0
  583. DO 433 L=1,NPG
  584. UH=0.D0
  585. DO 4333 KH=1,IDIM
  586. UH=UH+HR(KH,J,L)*HR(KH,I,L)
  587. 4333 CONTINUE
  588. U=U-UH*PGSQ(L)*DEUPI*BETA*IZTGG1.VPOCHA(KC,1)
  589. 433 CONTINUE
  590. AM(NK,I,J)=U
  591. 44 CONTINUE
  592.  
  593. 43 CONTINUE
  594.  
  595. ELSEIF(IK1.EQ.4.OR.IK1.EQ.5)THEN
  596. C Cas coefficient SOMMET
  597.  
  598. DO 45 KK=1,NBEL
  599. NK=NK+1
  600. DO 46 I=1,NP
  601. J=MELEMP.NUM(I,NK)
  602. DO 46 N=1,IDIM
  603. XYZ(N,I)=XCOOR((J-1)*(IDIM+1) +N)
  604. 46 CONTINUE
  605.  
  606. CALL CALJBR
  607. &(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
  608.  
  609. DO 47 I=1,NP
  610. DO 47 J=1,NP
  611. U=0.D0
  612. DO 453 L=1,NPG
  613. UH=0.D0
  614. DO 4533 KH=1,IDIM
  615. UH=UH+HR(KH,J,L)*HR(KH,I,L)
  616. 4533 CONTINUE
  617.  
  618. UA=0.D0
  619. DO 454 JJ=1,NP
  620. J1=LECT(IPT1.NUM(JJ,NK))
  621. DO 454 KH=1,IDIM
  622. KG=(IK1-4)*(KH-1)+1
  623. UA=UA+FN(JJ,L)*IZTGG1.VPOCHA(J1,KG)
  624. 454 CONTINUE
  625. UA=UA/KG
  626.  
  627. U=U-UH*PGSQ(L)*DEUPI*BETA*UA
  628. 453 CONTINUE
  629. AM(NK,I,J)=U
  630. 47 CONTINUE
  631.  
  632. 45 CONTINUE
  633. ENDIF
  634.  
  635. SEGDES IPT1
  636. 430 CONTINUE
  637.  
  638. C write(6,*)'C Fin Stabilisation de toutes sortes'
  639. ENDIF
  640. C Fin Stabilisation de toutes sortes
  641.  
  642. NBME=IDIM
  643. CALL KRIPAD(MELEMI,MLENTI1)
  644. SEGACT MELEMI
  645. NBSOUS=MELEMI.LISOUS(/1)
  646. IF(NBSOUS.EQ.0)NBSOUS=1
  647. SEGINI IMATRI
  648.  
  649. IF(IKAS.EQ.2)THEN
  650. KSPGD=MELEMS
  651. KSPGP=MELEMC
  652. IRIGEL(2,1)=MELEMI
  653. IRIGEL(1,1)=MELEMP
  654. ELSE
  655. KSPGP=MELEMS
  656. KSPGD=MELEMC
  657. IRIGEL(1,1)=MELEMI
  658. IRIGEL(2,1)=MELEMP
  659. ENDIF
  660. SEGACT MELEMP
  661.  
  662. C write(6,*)' ds kmac melemp=',IRIGEL(1,1)
  663. C write(6,*)' ds kmac melemd=',IRIGEL(2,1)
  664.  
  665. IRIGEL(4,1)=IMATRI
  666. IF(NOMUL.EQ.0)THEN
  667. IRIGEL(7,1)=2
  668. ELSE
  669. IF(IKAS.EQ.1)IRIGEL(7,1)=-3
  670. IF(IKAS.EQ.2)IRIGEL(7,1)=3
  671. IF(IKAS.EQ.3)IRIGEL(7,1)=4
  672. ENDIF
  673.  
  674. SEGACT MELEMP
  675.  
  676. K0=0
  677. DO 11 L=1,MAX(1,MELEMI.LISOUS(/1))
  678. IPT1=MELEMI
  679. IPT2=MELEMP
  680. IF(MELEMI.LISOUS(/1).NE.0)IPT1=MELEMI.LISOUS(L)
  681. IF(MELEMP.LISOUS(/1).NE.0)IPT2=MELEMP.LISOUS(L)
  682. SEGACT IPT1,IPT2
  683. NBEL=IPT1.NUM(/2)
  684.  
  685. IF(IKAS.EQ.2)THEN
  686. MP=IPT1.NUM(/1)
  687. C avt NP=MELEMP.NUM(/1)
  688. NP=IPT2.NUM(/1)
  689. ELSE
  690. NP=IPT1.NUM(/1)
  691. C avt MP=MELEMP.NUM(/1)
  692. MP=IPT2.NUM(/1)
  693. ENDIF
  694.  
  695. DO 12 I=1,NBME
  696. SEGINI IZAFM
  697. LIZAFM(L,I)=IZAFM
  698. IF(IKAS.EQ.2)THEN
  699. WRITE(NOM,FMT='(I1,A3)')I,NOMD(1:3)
  700. LISDUA(I)=NOM//' '
  701. LISPRI(I)=NOMP
  702. ELSE
  703. WRITE(NOM,FMT='(I1,A3)')I,NOMP(1:3)
  704. LISPRI(I)=NOM//' '
  705. LISDUA(I)=NOMD
  706. ENDIF
  707. 12 CONTINUE
  708. IPM1=LIZAFM(L,1)
  709. IPM2=LIZAFM(L,2)
  710. IPM3=LIZAFM(L,2)
  711. IF(IDIM.EQ.3)IPM3=LIZAFM(L,3)
  712.  
  713. IF(IK1.LT.4)THEN
  714. C write(6,*)' Appel KPRUSS'
  715. CALL KPRUSS(IPT1,IPT2,
  716. &IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,IK1,K0,IKOMP)
  717. C write(6,*)'Retour KPRUSS'
  718. ELSE
  719. C write(6,*)' Appel KPRJSS'
  720. CALL KPRJSS(IPT1,IPT2,IPM1,IPM2,IPM3,
  721. & IAXI,IKAS,INEFMD,KPRE,IZTGG1,MLENTI,IK1,IKOMP)
  722. ENDIF
  723.  
  724. K0=K0+NBEL
  725. SEGDES IPM1,IPM2,IPM3
  726. SEGDES IPT1
  727. 11 CONTINUE
  728. SEGDES MELEMI
  729. IF(IK1.EQ.4)SEGSUP MLENTI
  730.  
  731.  
  732. C%%%%%%%%%%%%%%%%% CAS DES PRESSIONS CONTINUES %%%%%%%%%%%%%%%%%%%%%%%
  733.  
  734. IF(KPRE.EQ.5)THEN
  735. CALL LEKTAB(MTABZ,'MMAIL',MELEME)
  736. CALL LEKTAB(MTABZ,'MSOMMET',MELEMS)
  737.  
  738. SEGACT MELEME
  739. NBSOU1=LISOUS(/1)
  740. IF(NBSOU1.EQ.0)NBSOU1=1
  741.  
  742. IRIGEL(1,2)=MELEME
  743. IRIGEL(2,2)=MELEME
  744. IRIGEL(7,2)=0
  745.  
  746. c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  747.  
  748. IHV=0
  749. AIMPL=1.D0
  750. C --- CALCUL DU NOMBRE DE PAQUETS DE LRV ELEMENTS ---
  751. CALL KPLRVM(MELEME,LRV,NBSOUS)
  752. NBME=1
  753.  
  754. SEGINI IMATRI
  755. IRIGEL(4,2)=IMATRI
  756. KSPGP=MELEMS
  757. KSPGD=MELEMS
  758. LISPRI(1)=NOMD//' '
  759. LISDUA(1)=NOMD//' '
  760.  
  761. LL=0
  762. NUTOEL=0
  763. SEGACT MELEME
  764.  
  765. DO 101 L=1,MAX(1,LISOUS(/1))
  766. IPT1=MELEME
  767. IF(LISOUS(/1).NE.0)IPT1=LISOUS(L)
  768. SEGACT IPT1
  769. NOM0=NOMS(IPT1.ITYPEL)//' '
  770. CALL KALPBG(NOM0,'FONFORM ',IZFFM)
  771. SEGACT IZFFM*MOD
  772. IZHR=KZHR(1)
  773. SEGACT IZHR*MOD
  774. NES=GR(/1)
  775. NPG=GR(/3)
  776.  
  777. NP = IPT1.NUM(/1)
  778. MP = NP
  779. NBELEM=IPT1.NUM(/2)
  780.  
  781. NNN=MOD(NBELEM,LRV)
  782. IF(NNN.EQ.0) NPACK=NBELEM/LRV
  783. IF(NNN.NE.0) NPACK=1+(NBELEM-NNN)/LRV
  784.  
  785. DO 701 KPACK=1,NPACK
  786.  
  787. C --- POUR CHAQUE PAQUET DE LRV ELEMENTS ou moins
  788. KDEB=1+(KPACK-1)*LRV
  789. KFIN=MIN(NBELEM,KDEB+LRV-1)
  790.  
  791. NBEL=KFIN-KDEB+1
  792. LL=LL+1
  793.  
  794. SEGINI IPM1,IPS1
  795. C write(6,*)' IPM1,LL=',IPM1,LL
  796. LIZAFM(LL,1)=IPM1
  797.  
  798. IPM2=IPM1
  799. IPM3=IPM1
  800. IPS2=IPS1
  801. IPS3=IPS1
  802.  
  803. KITT=2
  804. KJTT=IK1
  805. VISCO=IZTGG1
  806. CALL INITD(VISC,3, 1.D0)
  807.  
  808. CALL XLAPL(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,
  809. & VISC,VISC,VISC,KITT,KJTT,IK1,
  810. & IPT1.NUM(1,KDEB),NBEL,NUTOEL,XCOOR,AIMPL,IKOMP,
  811. & IPM1.AM,IPM2.AM,IPM3.AM,
  812. & IPS1.AM,IPS2.AM,IPS3.AM,
  813. & NINKO,IHV,IARG,VISC)
  814.  
  815. NUTOEL=NUTOEL+NBEL
  816. SEGDES IPM1
  817.  
  818. 701 CONTINUE
  819. SEGDES IPT1
  820. 101 CONTINUE
  821.  
  822.  
  823.  
  824.  
  825.  
  826.  
  827. c%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  828.  
  829. ENDIF
  830.  
  831. CALL ECROBJ('MATRIK',MATRIK)
  832.  
  833. NAT=2
  834. NSOUPO=0
  835. SEGINI MCHPOI
  836. JATTRI(1)=2
  837. CALL ECROBJ('CHPOINT',MCHPOI)
  838.  
  839. C write(6,*)' Fin operateur KMIC'
  840. SEGDES IMATRI,MATRIK
  841. RETURN
  842. 1001 FORMAT(20(1X,I5))
  843. 1002 FORMAT(10(1X,1PE11.4))
  844. END
  845.  
  846.  
  847.  
  848.  
  849.  
  850.  
  851.  
  852.  
  853.  
  854.  
  855.  
  856.  
  857.  
  858.  
  859.  
  860.  
  861.  

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