Télécharger yfrt.eso

Retour à la liste

Numérotation des lignes :

yfrt
  1. C YFRT SOURCE CB215821 24/04/12 21:17:32 11897
  2. SUBROUTINE YFRT(NOMPR,MTABX,IHV,MZIN,NOMI,MATRIK,MCHPO1)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. C***********************************************************************
  6. C
  7. C CET OPERATEUR DISCRETISE Convection Diffusion + source
  8. C
  9. C SYNTAXE :
  10. C ---------
  11. C
  12. C
  13. C COEFFICIENTS :
  14. C --------------
  15. C
  16. C
  17. C ALF (SCAL DOMA) DIFFUSIVITE THERMIQUE
  18. C (SCAL ELEM)
  19. C (SCAL NOEU)
  20. C
  21. C INCONNUES :
  22. C -----------
  23. C
  24. C TN CHAMP DE TEMPERATURE
  25. C UN CHAMP DE VITESSE TRANSPORTANT
  26. C
  27. C NOMPR Nom de l'opérateur appelant (NS,TSCA,KONV ou LAPN)
  28. C suivant l'opérateur les coefs sont différents
  29. C Schéma en temps ou l'opérateur intervient :
  30. C SEMI CNG TVISQ
  31. C teta - schéma Crank Nicolson généralisé Tenseur visqueux
  32. C Ces schéma nécessitent que DFDT soit en CENTREE
  33. C Le Schéma BDF2 ne fait intervenir que DFDT
  34. C
  35. C E/ MTABX Table de l'opérateur -> coefficients
  36. C E/ IHV=0 SCALAIRE
  37. C IHV=1 VITESSE
  38. C E/ MZIN pointeur du CHPOINT de l'inconnue
  39. C E/ NOMI nom de l'inconnue
  40. C************************************************************************
  41.  
  42.  
  43. -INC PPARAM
  44. -INC CCOPTIO
  45. -INC CCREEL
  46. -INC CCGEOME
  47. -INC SMCOORD
  48. -INC SMLENTI
  49. POINTEUR MLENT4.MLENTI
  50. -INC SMLMOTS
  51. -INC SMMODEL
  52. -INC SMELEME
  53. POINTEUR IGEOM.MELEME,MELEMS.MELEME,MELEMC.MELEME,MELEMP.MELEME
  54. POINTEUR MELEM2.MELEME
  55. -INC SMCHPOI
  56. POINTEUR MZIN.MCHPOI
  57. POINTEUR MTETA1.MPOVAL,MTETA2.MPOVAL,MTETA3.MPOVAL,MTETA4.MPOVAL
  58. -INC SMTABLE
  59. POINTEUR MTABZ.MTABLE
  60. -INC SIZFFB
  61. POINTEUR IZF1.IZFFM,IZH2.IZHR,IZW.IZFFM,IZWH.IZHR
  62. SEGMENT SAJT
  63. REAL*8 AJT(IDIM,IDIM,NPG),RF1(NP,MP,IDIM),SM1(NP,IDIM)
  64. REAL*8 TN1(NP,IDIM),TN2(NP,IDIM)
  65. REAL*8 WT(MP,NPG), WS(MP,NPG)
  66. REAL*8 WTC(IDIM*NPG,MP),WSC(IDIM*NPG,MP),HRD(IDIM*NPG,NP)
  67. REAL*8 FND(NPG,NP)
  68. ENDSEGMENT
  69. -INC SMRIGID
  70. C-INC SMMATRIK
  71. C*******************************************************************
  72. C
  73. SEGMENT MATRIK
  74. REAL*8 COEMTK(NMATRI)
  75. INTEGER JRIGEL(NRIGE,NMATRI)
  76. INTEGER KSYM,KMINC,KMINCP,KMINCD,KIZM
  77. INTEGER KISPGT,KISPGP,KISPGD
  78. INTEGER KNTTT,KNTTP,KNTTD
  79. INTEGER KIDMAT(NKID)
  80. INTEGER KKMMT(NKMT)
  81. ENDSEGMENT
  82.  
  83. SEGMENT JMATRI
  84. CHARACTER*8 LISPRI(NBMF),LJSDUA(NBMF)
  85. INTEGER LIZAFM(NBSOUS,NBMF)
  86. INTEGER KSPGP,KSPGD
  87. ENDSEGMENT
  88. POINTEUR JMATRS.JMATRI,JMATR1.JMATRI,JMATR2.JMATRI,JMATR3.JMATRI
  89.  
  90. C Stokage matrices elementaires non assemblees (valeurs)
  91. SEGMENT IZAFM
  92. REAL*8 AM(NBEL,NP,MP)
  93. ENDSEGMENT
  94. POINTEUR IPM1.IZAFM,IPM2.IZAFM,IPM3.IZAFM,IPM4.IZAFM
  95.  
  96. C*******************************************************************
  97. POINTEUR IPM.IZAFM
  98. -INC SMCHAML
  99.  
  100. LOGICAL TDFDT,TCONV,TSOUR,TSOURB,TMATRI
  101. LOGICAL ICAL,XPG,XRIG
  102. LOGICAL XDIAG,XTV,XTG,XBDF,XCONS
  103. CHARACTER*(*) NOMPR,NOMI
  104. CHARACTER*8 CHHH,TYPE,NOM,NOM0,CHAI,SCHT,NOMPER,MPRE
  105. CHARACTER*4 INCOD
  106. CHARACTER*(LOCOMP) NOMP(3),NOMD(3),NMACO
  107. DIMENSION ICOEF(4),SU(3),XPOI(3)
  108. C*****************************************************************************
  109. CTCLSF
  110. NOMPER=NOMPR
  111. write(6,*)' DEBUT YFRT appele par ',NOMPER
  112.  
  113. C- Table de l'opérateur (pointeur MTABX)
  114.  
  115. C- Récupération de la table RV (pointeur MTAB1)
  116. CALL LEKTAB(MTABX,'EQEX',MTAB1)
  117. C- Récupération de la table des options de l'opérateur (pointeur KOPTI)
  118. CALL LEKTAB(MTABX,'KOPT',KOPTI)
  119. C- Récupération de la table DOMAINE de la zone (pointeur MTABZ)
  120. CALL LEKTAB(MTABX,'DOMZ',MTABZ)
  121. C- Récupération de la table INCO (pointeur KINC)
  122. CALL LEKTAB(MTAB1,'INCO',KINC)
  123. IF(IERR.NE.0) RETURN
  124.  
  125. C*****************************************************************************
  126. C Traitement des options
  127. C KIMPL = 0 -> EXPL 1 -> IMPL 2 -> CN
  128. C AIMPL=0 explicite
  129. C KKALE=1 Indice indiquant que l'on est en ALE
  130. C sinon KKALE=0
  131. C KFORM = 0 -> SI 1 -> EF 2 -> VF 3 -> EFMC
  132. C IDCEN : entier indiquant le type de décentrement souhaité
  133. C IDCEN 1-> CENTREE 2-> SUPGDC 3-> SUPG 4-> TVISQUEU 5-> CNG
  134. C E/ CMD : Coefficient multiplicateur du décentrement
  135. C Si IDCEN=4 ou =5 CMD=DT
  136.  
  137. CALL ACME(KOPTI,'KFORM',KFORM)
  138. CALL ACME(KOPTI,'KIMPL',KIMPL)
  139. CALL ACMF(KOPTI,'AIMPL',AIMPL)
  140. IKOMP=0
  141. *? CALL ACME(KOPTI,'IKOMP',IKOMP)
  142. CALL ACME(KOPTI,'KMACO',KMACO)
  143. KKMACO=KMACO
  144. CALL ACME(KOPTI,'IDCEN',IDCEN)
  145. CALL ACMF(KOPTI,'CMD',CMD)
  146. *? CALL ACME(KOPTI,'ALE',KKALE)
  147.  
  148. C Restrictions des options
  149. IF(NOMPER.NE.'NS')KKALE=0
  150. IF(NOMPER.EQ.'LAPN')IDCEN=1
  151.  
  152. IF(KIMPL.EQ.0)AIMPL=0.D0
  153. IF(KIMPL.EQ.1)AIMPL=1.D0
  154.  
  155. c write(6,*)' YFRT KMACO=',KMACO,'IDCEN=',IDCEN
  156. IF(KMACO.EQ.1.AND.IDCEN.NE.5)THEN
  157. CALL ACMM(KOPTI,'NMACO',NMACO)
  158. TYPE=' '
  159. CALL ACMO(MTAB1,NMACO,TYPE,MATRIK)
  160. IF(TYPE.EQ.'MATRIK')THEN
  161. KKMACO=2
  162. ELSE
  163. KKMACO=1
  164. ENDIF
  165. ENDIF
  166. IF(IDCEN.EQ.5)KKMACO=0
  167. c write(6,*)' YFRT KMACO=',kmaco,' KKMACO=',KKMACO,'IKOMP=',ikomp
  168.  
  169. C*****************************************************************************
  170.  
  171. KPOIND=0
  172. NBMF=1
  173. XRIG=.FALSE.
  174. TMATRI=.TRUE.
  175. XPG=.FALSE.
  176. IF(IDCEN.GT.1)XPG=.TRUE.
  177. IF(IDCEN.EQ.4)AIMPL=0.D0
  178. IF(IDCEN.EQ.5)THEN
  179. AIMPL=1.D0
  180. ENDIF
  181.  
  182. AIEX=AIMPL
  183. IF(AIMPL.EQ.0.D0.AND.KKMACO.NE.0)THEN
  184. TMATRI=.TRUE.
  185. AIEX=1.D0
  186. ENDIF
  187.  
  188. IF(AIMPL.EQ.0.D0.AND.KKMACO.EQ.0)THEN
  189. TMATRI=.FALSE.
  190. ENDIF
  191.  
  192. IAXI=0
  193. IF(IFOMOD.EQ.0)IAXI=2
  194. DEUPI=1.D0
  195. IF(IAXI.NE.0)DEUPI=2.D0*XPI
  196.  
  197. IF(XRIG)THEN
  198. IF(IHV.EQ.0)THEN
  199. NOMP(1)='T '
  200. NOMD(1)='Q '
  201. ELSEIF(IHV.EQ.1)THEN
  202. NOMP(1)='UX '
  203. NOMP(2)='UY '
  204. NOMP(3)='UZ '
  205. NOMD(1)='FX '
  206. NOMD(2)='FY '
  207. NOMD(3)='FZ '
  208. ENDIF
  209. ELSE
  210. IF(IHV.EQ.0)THEN
  211. NOMP(1)=NOMI
  212. NOMD(1)=NOMI
  213. ELSEIF(IHV.EQ.1)THEN
  214. WRITE(NOMP(1),FMT='(I1,A3)')1,NOMI(1:LOCOMP-1)
  215. WRITE(NOMD(1),FMT='(I1,A3)')1,NOMI(1:LOCOMP-1)
  216. WRITE(NOMP(2),FMT='(I1,A3)')2,NOMI(1:LOCOMP-1)
  217. WRITE(NOMD(2),FMT='(I1,A3)')2,NOMI(1:LOCOMP-1)
  218. WRITE(NOMP(3),FMT='(I1,A3)')3,NOMI(1:LOCOMP-1)
  219. WRITE(NOMD(3),FMT='(I1,A3)')3,NOMI(1:LOCOMP-1)
  220. ENDIF
  221. ENDIF
  222. c write(6,*)' NOMP=',nomp
  223. c write(6,*)' NOMD=',nomd
  224.  
  225. CALL LEKTAB(MTABZ,'SOMMET',MELEMS)
  226. CALL LEKTAB(MTABZ,'MAILLAGE',MELEME)
  227. MELEM2=MELEME
  228.  
  229. C On cree un second membre non vide
  230. SEGACT MELEMS
  231. N=MELEMS.NUM(/2)
  232. IF(IHV.EQ.0)NC=1
  233. IF(IHV.EQ.1)NC=IDIM
  234. SEGDES MELEMS
  235. NSOUPO=1
  236. NAT=2
  237. SEGINI MCHPO1,MSOUP1,MPOVA1
  238. MCHPO1.JATTRI(1)=2
  239. MCHPO1.IFOPOI=IFOUR
  240. MCHPO1.MTYPOI=' '
  241. MCHPO1.MOCHDE=' '
  242. MCHPO1.IPCHP(1)=MSOUP1
  243. MSOUP1.IGEOC=MELEMS
  244. MSOUP1.IPOVAL=MPOVA1
  245. IF(IHV.EQ.0)THEN
  246. MSOUP1.NOCOMP(1)=NOMD(1)
  247. ELSE
  248. DO 91 N=1,NC
  249. MSOUP1.NOCOMP(N)=NOMD(N)
  250. 91 CONTINUE
  251. ENDIF
  252. SEGDES MSOUP1,MCHPO1
  253.  
  254. C*************************************************************************
  255. C Lecture des Inco(s) aux temps precedents si Transitoire
  256. c write(6,*)'Lecture des INCO MZIN=',MZIN
  257.  
  258. CALL LRCHT(MZIN,MTETA1,TYPE,IGEOM)
  259. CALL KRIPAD(IGEOM,MLENT1)
  260. CALL VERPAD(MLENT1,MELEMS,IRET)
  261. IF(IRET.NE.0)THEN
  262. C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon support géométrique
  263. MOTERR(1: 8) = 'TETA1'
  264. MOTERR(9:16) = 'CHPOINT '
  265. WRITE(IOIMP,*)'Opérateur : ',NOMPER
  266. CALL ERREUR(788)
  267. RETURN
  268. ENDIF
  269.  
  270. NC=MTETA1.VPOCHA(/2)
  271.  
  272. IF((IHV.EQ.0.AND.NC.NE.1).OR.
  273. & (IHV.EQ.1.AND.NC.NE.IDIM))THEN
  274. C% Indice %m1:8 : L'objet %m9:16 n'a pas le bon nombre de composantes
  275. MOTERR( 1: 8) = 'Inconnue'
  276. MOTERR( 9:16) = 'CHPOINT'
  277. WRITE(IOIMP,*)'Operateur : ',NOMPER
  278. CALL ERREUR(784)
  279. RETURN
  280. ENDIF
  281.  
  282. SEGDES IGEOM,MELEMS
  283.  
  284. C*****************************************************************************
  285. C Lecture du coefficient
  286. c write(6,*)'Lecture des COEF '
  287.  
  288. CALL ACME(MTABX,'IARG',IARG)
  289. c write(6,*)' YFRT IARG=',IARG,' KKALE=',KKALE
  290.  
  291. C--------- CAS FROT Formulation non conservative ----------------------------
  292. IF(NOMPER.EQ.'FROT')THEN
  293.  
  294. C Coef K={Kx Ky Kz}
  295. NUCOEF=1
  296. c write(6,*)' Coef 1 ',NUCOEF,' MTABZ=',MTABZ
  297. CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,1,MCHEL1,KPOIND,0,MCHELG)
  298. IF (IERR.NE.0) RETURN
  299.  
  300. C Coef Beta={Betax Betay Betaz}
  301. NUCOEF=NUCOEF+1
  302. c write(6,*)' Coef 2 ',NUCOEF,' MTABZ=',MTABZ
  303. CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,1,MCHEL2,KPOIND,0,MCHELG)
  304. IF (IERR.NE.0) RETURN
  305.  
  306. C Un
  307. NUCOEF=NUCOEF+1
  308. c write(6,*)' Coef 3 ',NUCOEF
  309. CALL LEKMOF(MTABZ,NUCOEF,MTABX,KINC,1,MCHEL3,KPOIND,0,MCHELG)
  310. IF (IERR.NE.0) RETURN
  311.  
  312.  
  313. c write(6,*)'FIN CAS FROT Formulation non conservative'
  314. C----- FIN CAS FROT Formulation non conservative ----------------------------
  315.  
  316.  
  317. ENDIF
  318. c write(6,*)' Fin lect COEF '
  319. C*************************************************************************
  320. C******* CALCUL **********************************************************
  321. C*************************************************************************
  322.  
  323. SEGACT MTETA1
  324.  
  325. CALL KRIPAD(MELEMS,MLENTI)
  326. SEGDES MELEMS
  327.  
  328. SEGACT MELEME,MELEM2
  329. NBSOUS=LISOUS(/1)
  330. IF(NBSOUS.EQ.0)NBSOUS=1
  331.  
  332. C MATRIK . MATRIK . MATRIK . MATRIK . MATRIK . MATRIK MATRIK . MATRIK .
  333. IF(KKMACO.NE.2)THEN
  334. C ... Création / Ajustement des matrices élémentaires
  335. IF(TMATRI)THEN
  336. C Cas RIGIDITE
  337. IF(XRIG)THEN
  338.  
  339. NRIGE=8
  340. NRIGEL=0
  341. SEGINI MRIGID
  342.  
  343. NRIGE=8
  344. NRI=IRIGEL(/2)
  345. NRIGEL=NBSOUS+NRI
  346. SEGADJ MRIGID
  347. c write(6,*)' NRIGE,NRIGEL,MRIGID=',NRIGE,NRIGEL,MRIGID
  348.  
  349. ELSE
  350. C Cas MATRIK
  351.  
  352. NRIGE=7
  353. NKID =9
  354. NKMT =7
  355. NMATRI=0
  356. SEGINI MATRIK
  357. c write(6,*)' YFRT creation MATRIK =',matrik
  358.  
  359. NRIGE=7
  360. NKID =9
  361. NKMT =7
  362. NMATR0=JRIGEL(/2)
  363.  
  364. NBME=1
  365. IF(IHV.EQ.1)NBME=IDIM
  366. NMATRI=NBME
  367. SEGADJ MATRIK
  368.  
  369. DO 41 M=1,NBME
  370. JRIGEL(1,M)=MELEME
  371. JRIGEL(2,M)=MELEM2
  372. JRIGEL(7,M)=0
  373. c? IF(TCONV)JRIGEL(7,M)=2
  374. SEGINI JMATRI
  375. JRIGEL(4,M)=JMATRI
  376. KSPGP=MELEMS
  377. KSPGD=MELEMS
  378. LISPRI(1)=NOMP(M)//' '
  379. LJSDUA(1)=NOMD(M)//' '
  380. 41 CONTINUE
  381. ENDIF
  382. ENDIF
  383. C ... Fin Création / Ajustement des matrices élémentaires
  384.  
  385.  
  386. ELSE
  387. C Activation RIGIDITE
  388. IF(XRIG)THEN
  389. SEGACT MRIGID
  390. ELSE
  391. C Activation MATRIK
  392. SEGACT MATRIK
  393. ENDIF
  394.  
  395. C ... Activation des matrices élémentaires
  396. IF(TMATRI)THEN
  397. C Cas RIGIDITE
  398. IF(XRIG)THEN
  399. NRIGE=8
  400. NRI=IRIGEL(/2)
  401. NRIGEL=NBSOUS+NRI
  402. SEGAct MRIGID
  403. c write(6,*)' NRIGE,NRIGEL,MRIGID=',NRIGE,NRIGEL,MRIGID
  404.  
  405. ELSE
  406. C Cas MATRIK Pleine
  407. NMATRI=JRIGEL(/2)
  408. NBME=NMATRI
  409.  
  410. DO 42 M=1,NBME
  411. JMATRI=JRIGEL(4,M)
  412. SEGACT JMATRI*MOD
  413. LISPRI(1)=NOMP(M)//' '
  414. LJSDUA(1)=NOMD(M)//' '
  415. 42 CONTINUE
  416.  
  417. ENDIF
  418. ENDIF
  419. C ... Fin Activation des matrices élémentaires
  420.  
  421. ENDIF
  422. C FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN M
  423. C ......................................................................
  424.  
  425. SEGACT MCHEL1
  426. SEGACT MCHEL2
  427. SEGACT MCHEL3
  428.  
  429. IF(MAX(1,MELEM2.LISOUS(/1)).NE.MAX(1,LISOUS(/1)))THEN
  430. WRITE(6,*)' Geometries incompatibles dans ',nomper
  431. C% Données incompatibles
  432. CALL ERREUR(22)
  433. RETURN
  434. ENDIF
  435.  
  436. NKD=0
  437. DO 101 L=1,MAX(1,LISOUS(/1))
  438. IPT1=MELEME
  439. IPT2=MELEM2
  440. IF(LISOUS(/1).NE.0)IPT1=LISOUS(L)
  441. IF(MELEM2.LISOUS(/1).NE.0)THEN
  442. IPT2=MELEM2.LISOUS(L)
  443. NKD=0
  444. ENDIF
  445. SEGACT IPT1,IPT2
  446.  
  447. NOM0 = NOMS(IPT1.ITYPEL)//' '
  448. CALL KALPBG(NOM0,'FONFORM ',IZFFM)
  449. IF(IZFFM.EQ.0)RETURN
  450. SEGACT IZFFM*MOD
  451. IZHR=KZHR(1)
  452. SEGACT IZHR*MOD
  453. IZF1 = KTP(1)
  454. IZH2 = KZHR(2)
  455. IZW = IZFFM
  456. IZWH = IZHR
  457.  
  458. NES=GR(/1)
  459. NPG=GR(/3)
  460.  
  461. NP = IPT1.NUM(/1)
  462. MP = IPT2.NUM(/1)
  463. C? MP = IZW.FN(/1) ceci doit etre identique
  464.  
  465. NBEL=IPT1.NUM(/2)
  466.  
  467. SEGINI SAJT
  468.  
  469. c.......................................................................
  470. C MATRIK . MATRIK . MATRIK . MATRIK . MATRIK . MATRIK MATRIK . MATRIK .
  471. c write(6,*)' KKMACO=',KKMACO,'TMATRI=',TMATRI
  472. IF(KKMACO.NE.2)THEN
  473. IF(TMATRI)THEN
  474. C Cas RIGIDITE
  475. IF(XRIG)THEN
  476. IRIGEL(1,NRI+L)=MELEME
  477. COERIG(L)=1.D0
  478.  
  479. IRIGEL(7,NRI+L)=0
  480. c? IF(TCONV)IRIGEL(7,NRI+L)=2
  481.  
  482. NBME=1
  483. IF(IHV.EQ.1)NBME=IDIM
  484. NLIGRP=NP
  485. NLIGRD=MP
  486. SEGINI DESCR
  487. IRIGEL(3,NRI+L)=DESCR
  488. IF(NBME.EQ.1)THEN
  489. DO 102 I=1,NLIGRP
  490. LISINC(I)=NOMP(1)//' '
  491. NOELEP(I)=I
  492. 102 CONTINUE
  493. DO 103 I=1,NLIGRD
  494. LISDUA(I)=NOMD(1)//' '
  495. NOELED(I)=I
  496. 103 CONTINUE
  497. ELSE
  498. ENDIF
  499. SEGDES DESCR
  500.  
  501. NELRIG=NBEL
  502. SEGINI xMATRI
  503. IRIGEL(4,NRI+L)=xMATRI
  504. xmatri.symre=irigel(7,nri+l)
  505. c write(6,*)'NELRIG,IMATRI=',NELRIG,IMATRI
  506.  
  507. * DO 104 K=1,NELRIG
  508. * SEGINI XMATRI
  509. c write(6,*)'NLIGRD,NLIGRP,XMATRI=',NLIGRD,NLIGRP,XMATRI
  510. * IMATTT(K)=XMATRI
  511. * 104 CONTINUE
  512.  
  513. ELSE
  514. C Cas MATRIK
  515. SEGINI IPM1
  516. c mtail=(IPM1.AM(/1))*(ipm1.am(/2))*(ipm1.am(/3))
  517. JMATR1=JRIGEL(4,NMATR0+1)
  518. JMATR1.LIZAFM(L,1)=IPM1
  519. IPM2=IPM1
  520. IPM3=IPM1
  521. IF(NBME.GE.2)THEN
  522. JMATR2=JRIGEL(4,NMATR0+2)
  523. SEGINI IPM2
  524. c mtail=(IPM2.AM(/1))*(ipm2.am(/2))*(ipm2.am(/3))
  525. JMATR2.LIZAFM(L,1)=IPM2
  526. ICAL=.TRUE.
  527. ENDIF
  528. IF(NBME.GE.3)THEN
  529. JMATR3=JRIGEL(4,NMATR0+3)
  530. SEGINI IPM3
  531. JMATR3.LIZAFM(L,1)=IPM3
  532. ICAL=.TRUE.
  533. ENDIF
  534. ENDIF
  535. ENDIF
  536. c write(6,*)' YFRT appelé par ',NOMPER,' Taille matrice=',mtail
  537.  
  538. ELSE
  539.  
  540. C Cas RIGIDITE
  541. IF(XRIG)THEN
  542.  
  543. ELSE
  544. C Cas MATRIK
  545. JMATR1=JRIGEL(4,1)
  546. SEGACT JMATR1
  547. IPM1=JMATR1.LIZAFM(L,1)
  548. IPM2=IPM1
  549. IPM3=IPM1
  550. SEGACT IPM1
  551. IF(NBME.GE.2)THEN
  552. JMATR2=JRIGEL(4,2)
  553. SEGACT JMATR2
  554. IPM2=JMATR2.LIZAFM(L,1)
  555. SEGACT IPM2
  556. ENDIF
  557. IF(NBME.GE.3)THEN
  558. JMATR3=JRIGEL(4,3)
  559. IPM3=JMATR3.LIZAFM(L,1)
  560. SEGACT IPM3
  561. ENDIF
  562. ENDIF
  563. ENDIF
  564. C FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN MATRIK . FIN M
  565. C ......................................................................
  566.  
  567.  
  568. C----K
  569. IK1=1
  570. MCHAM1=MCHEL1.ICHAML(L)
  571. SEGACT MCHAM1
  572.  
  573. MELVA1=MCHAM1.IELVAL(1)
  574. SEGACT MELVA1
  575. N1PTEL=MELVA1.VELCHE(/1)
  576. N1EL=MELVA1.VELCHE(/2)
  577. IF(N1EL.EQ.1)THEN
  578. IK1=1
  579. ELSEIF(N1EL.EQ.NBEL)THEN
  580. IK1=0
  581. ENDIF
  582.  
  583. C----Beta
  584. IK2=1
  585. MCHAM2=MCHEL2.ICHAML(L)
  586. SEGACT MCHAM2
  587.  
  588. MELVA2=MCHAM2.IELVAL(1)
  589. SEGACT MELVA2
  590. N1PTEL=MELVA2.VELCHE(/1)
  591. N1EL=MELVA2.VELCHE(/2)
  592. IF(N1EL.EQ.1)THEN
  593. IK2=1
  594. ELSEIF(N1EL.EQ.NBEL)THEN
  595. IK2=0
  596. ENDIF
  597.  
  598. C----U
  599. IK3=1
  600. MCHAM3=MCHEL3.ICHAML(L)
  601. SEGACT MCHAM3
  602.  
  603. MELVA3=MCHAM3.IELVAL(1)
  604. SEGACT MELVA3
  605. N1PTEL=MELVA3.VELCHE(/1)
  606. N1EL=MELVA3.VELCHE(/2)
  607. IF(N1EL.EQ.1)THEN
  608. IK3=1
  609. ELSEIF(N1EL.EQ.NBEL)THEN
  610. IK3=0
  611. ENDIF
  612.  
  613.  
  614. c write(6,*)' AVANT 108 NC=',NC,' NBEL=',NBEL,MP,NP,NC
  615. C===============================================
  616. AI1=AIMPL-1.D0
  617. IF(AIMPL.EQ.0.D0)THEN
  618. AI2=-1.D0
  619. ELSE
  620. AI2=AI1/AIMPL
  621. ENDIF
  622.  
  623. DO 108 KE=1,NBEL
  624.  
  625. NK1=KE + IK1*(1 - KE)
  626. NK2=KE + IK2*(1 - KE)
  627. NK3=KE + IK3*(1 - KE)
  628.  
  629. DO I=1,NP
  630. J=IPT1.NUM(I,KE)
  631. DO N=1,IDIM
  632. XYZ(N,I)=XCOOR((J-1)*(IDIM+1)+N)
  633. ENDDO
  634. ENDDO
  635.  
  636. DO I=1,NP
  637. I1=MLENT1.LECT(IPT1.NUM(I,KE))
  638. DO N=1,NC
  639. TN1(I,N)=MTETA1.VPOCHA(I1,N)
  640. ENDDO
  641. ENDDO
  642.  
  643. CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,
  644. * IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
  645.  
  646. IF (IDCEN.EQ.0.OR.IDCEN.EQ.1) THEN
  647. CALL RSETD(WT,FN,(NP*NPG))
  648. CALL RSETD(WS,FN,(NP*NPG))
  649. ELSE
  650. CALL CALDEC(WT,WS,XYZ,GR,HR,FN,NES,IDIM,NP,NPG,AJT,
  651. & IDCEN,CMD,MELVA1.VELCHE(1,NK1),MELVA2.VELCHE(1,NK2),
  652. & MELVA3.VELCHE(1,NK3),TN1,NC,IKOMP,XREF,Aire,KE)
  653. ENDIF
  654.  
  655. CALL INITD(SM1,(MP*IDIM),0.D0)
  656.  
  657. C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  658. IF(KKMACO.EQ.2)THEN
  659.  
  660. C ...... Chargement Rigidite ou Matrik
  661. C Cas RIGIDITE
  662. IF(XRIG)THEN
  663. * XMATRI=IMATTT(KE)
  664. DO I=1,MP
  665. DO J=1,NP
  666. RE(I,J,ke)=RF1(J,I,1)
  667. ENDDO
  668. ENDDO
  669. * SEGDES XMATRI
  670. ELSE
  671. C Cas MATRIK
  672. DO N=1,NBME
  673. JMATR1=JRIGEL(4,N)
  674. IPM4=JMATR1.LIZAFM(L,1)
  675. DO I=1,MP
  676. DO J=1,NP
  677. RF1(J,I,N)=IPM4.AM(KE,J,I)
  678. ENDDO
  679. ENDDO
  680. ENDDO
  681. ENDIF
  682. C...... Multiplication Matrice Vecteur
  683.  
  684. DO 710 I=1,MP
  685. DO 711 J=1,NP
  686.  
  687. DO 716 N=1,NC
  688. SM1(I,N) = SM1(I,N) + AI2*RF1(J,I,N)*TN1(J,N)
  689. 716 CONTINUE
  690.  
  691. 711 CONTINUE
  692. 710 CONTINUE
  693.  
  694. ELSEIF(KKMACO.NE.2)THEN
  695. C..............................................................
  696. C...... Frot .... Formulation non conservative
  697. C préparation
  698.  
  699. DO 416 N=1,NC
  700. DO 410 I=1,MP
  701. c DO 411 J=1,NP
  702. U3=0.D0
  703. DO 412 LG=1,NPG
  704. PDR=PGSQ(LG)*DEUPI*RPG(LG)
  705. C1=MELVA1.VELCHE((N-1)*NPG+LG,NK1)*PDR
  706. C2=MELVA2.VELCHE((N-1)*NPG+LG,NK2)
  707. C3=ABS ( MELVA3.VELCHE((N-1)*NPG+LG,NK3) ) + 1.e-30
  708. CK=C1*(C3**(C2-1.D0))
  709. U3=U3+CK * FN(I,LG)
  710. 412 CONTINUE
  711.  
  712. SM1(I,N) = SM1(I,N) + AI1*U3*TN1(I,N)
  713. RF1(I,I,N) = AIEX*U3
  714. c? SM1(I,N) = SM1(I,N) + AI1*U3*TN1(J,N)
  715. c? RF1(J,I,N) = AIEX*U3
  716.  
  717. 411 CONTINUE
  718. 410 CONTINUE
  719. 416 CONTINUE
  720.  
  721.  
  722. C..............................................................
  723.  
  724.  
  725. ENDIF
  726.  
  727. C:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
  728. C=======================================================================
  729.  
  730. C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
  731. IF(TMATRI.AND.KKMACO.NE.2)THEN
  732. C ...... Chargement Rigidite ou Matrik
  733. C Cas RIGIDITE
  734. IF(XRIG)THEN
  735. * XMATRI=IMATTT(KE)
  736. DO I=1,NP
  737. DO J=1,NP
  738. RE(I,J,ke)=RF1(J,I,1)
  739. ENDDO
  740. ENDDO
  741. * SEGDES XMATRI
  742. ELSE
  743. C Cas MATRIK
  744. NBMM=1
  745. IF(ICAL)NBMM=NBME
  746. DO N=1,NBMM
  747. JMATR1=JRIGEL(4,NMATR0+N)
  748. IPM4=JMATR1.LIZAFM(L,1)
  749. DO I=1,NP
  750. DO J=1,NP
  751. IPM4.AM(KE,J,I)=RF1(J,I,N)
  752. ENDDO
  753. ENDDO
  754. ENDDO
  755. c write(6,*)' RF1 '
  756. c write(6,1002)(((RF1(J,I,N),j=1,np),i=1,np),n=1,nbmm)
  757. ENDIF
  758. ENDIF
  759.  
  760. C ...... Chargement Second membre
  761. DO I=1,NP
  762. I1=LECT(IPT1.NUM(I,KE))
  763. DO N=1,NC
  764. MPOVA1.VPOCHA(I1,N)=MPOVA1.VPOCHA(I1,N)+SM1(I,N)
  765. ENDDO
  766. ENDDO
  767. C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  768.  
  769. 108 CONTINUE
  770.  
  771. SEGSUP MCHAM1,MELVA1
  772. SEGSUP MCHAM2,MELVA2
  773. SEGSUP MCHAM3,MELVA3
  774.  
  775.  
  776.  
  777. SEGDES IPT1,IPT2
  778.  
  779. IF(TMATRI)THEN
  780. C Cas RIGIDITE
  781. IF(XRIG)THEN
  782. SEGDES xMATRI
  783. ELSE
  784. C Cas MATRIK
  785. SEGDES IPM1
  786. IF(NBME.GE.2)SEGDES IPM2
  787. IF(NBME.GE.3)SEGDES IPM3
  788. ENDIF
  789. ENDIF
  790.  
  791. SEGSUP IZFFM,IZHR,IZF1,IZH2
  792. SEGSUP SAJT
  793.  
  794. 101 CONTINUE
  795.  
  796. IF(TMATRI)THEN
  797. IF(.NOT.XRIG)THEN
  798. NBMM=JRIGEL(/2)
  799. IF(NBMM.NE.0)THEN
  800. DO 141 M=1,NBMM
  801. JMATRI=JRIGEL(4,M)
  802. SEGDES JMATRI
  803. 141 CONTINUE
  804. ENDIF
  805. ENDIF
  806. ENDIF
  807.  
  808. C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
  809. C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
  810. C"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
  811.  
  812.  
  813. SEGSUP MCHEL1
  814. SEGSUP MCHEL2
  815. SEGSUP MCHEL3
  816.  
  817. SEGDES MCHPO1,MPOVA1
  818. SEGDES MELEME,MELEM2
  819.  
  820. SEGSUP MLENTI
  821. c IF(TLINCO)THEN
  822. SEGSUP MLENT1
  823. SEGDES MTETA1
  824. c ENDIF
  825.  
  826. IF(TMATRI)THEN
  827. C Cas RIGIDITE
  828. IF(XRIG)THEN
  829. SEGDES MRIGID
  830. ELSE
  831. C Cas MATRIK
  832. SEGDES MATRIK
  833. ENDIF
  834. ENDIF
  835.  
  836.  
  837. c segact MCHPO1
  838. c MSOUP1=MCHPO1.IPCHP(1)
  839. c segact MSOUP1
  840. c write(6,*)MCHPO1.IPCHP(/1),MSOUP1.NOCOMP(/2)
  841. c write(6,*)(MSOUP1.NOCOMP(i),i=1,MSOUP1.NOCOMP(/2))
  842. c if(ihv.eq.1)call ecmo(mtabz,'TOTO','CHPOINT',MCHPO1)
  843.  
  844. C*************************************************************************
  845.  
  846. IF(KKMACO.EQ.1)THEN
  847. TYPE='MATRIK'
  848. CALL ECMO(MTAB1,NMACO,TYPE,MATRIK)
  849. ENDIF
  850.  
  851. IF(AIMPL.EQ.0.D0)THEN
  852. NRIGE=7
  853. NKID =9
  854. NKMT =7
  855. NMATRI=0
  856. SEGINI MATRIK
  857. ENDIF
  858.  
  859. c write(6,*)' FIN YFRT'
  860. RETURN
  861. 1001 FORMAT(20(1X,I5))
  862. 1002 FORMAT(10(1X,1PE11.4))
  863. END
  864.  
  865.  
  866.  
  867.  
  868.  
  869.  
  870.  
  871.  
  872.  
  873.  
  874.  
  875.  
  876.  
  877.  
  878.  
  879.  
  880.  
  881.  
  882.  
  883.  
  884.  
  885.  
  886.  
  887.  
  888.  
  889.  
  890.  
  891.  
  892.  
  893.  
  894.  
  895.  
  896.  
  897.  

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