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

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