Télécharger d2vtra.eso

Retour à la liste

Numérotation des lignes :

  1. C D2VTRA SOURCE BP208322 18/12/20 21:15:17 10048
  2. C
  3. SUBROUTINE D2VTRA(ITBAS,ITKM,ITA,KTKAM,IPMAIL,KTRES,KTNUM,KPREF,
  4. & KTPHI,KTLIAB,RIGIDE,ITCARA,LMODYN,IALGO,BETA)
  5. IMPLICIT INTEGER(I-N)
  6. IMPLICIT REAL*8(A-H,O-Z)
  7. *--------------------------------------------------------------------*
  8. * *
  9. * Operateur DYNE : algorithme de Fu - de Vogelaere *
  10. * ________________________________________________ *
  11. * *
  12. * Transpose l'information des objets de Castem2000 dans des *
  13. * tableaux de travail. *
  14. * *
  15. * Parametres: *
  16. * *
  17. * e ITBAS Table representant une base modale *
  18. * e ITKM Table contenant les matrices XK et XM *
  19. * e ITA Table contenant la matrice XASM *
  20. * es KTKAM Segment contenant les matrices XK, XASM et XM *
  21. * e IPMAIL Maillage de reference pour les CHPOINTs resultats *
  22. * es KTRES Segment de sauvegarde des resultats *
  23. * e KPREF Segment des points de reference *
  24. * es KTPHI Segment des deformees modales *
  25. * e KTLIAB Segment des liaisons sur base B *
  26. * e RIGIDE Vrai si corps rigide, faux sinon *
  27. * e IALGO =1 SI DE_VOGELAERE, =2 SI DIFFERENCES_CENTREES *
  28. * =3 si ACCELERATION_MOYENNE, =4 si FOX_GOODWIN *
  29. * *
  30. * Auteur, date de creation: *
  31. * *
  32. * Denis ROBERT-MOUGIN, le 26 mai 1989. *
  33. * *
  34. *--------------------------------------------------------------------*
  35. -INC SMRIGID
  36. -INC SMELEME
  37.  
  38. -INC PPARAM
  39. -INC CCOPTIO
  40. -INC CCREEL
  41. *
  42. SEGMENT,MTKAM
  43. REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
  44. REAL*8 XOPER(NB1,NB1,NOPER)
  45. ENDSEGMENT
  46. SEGMENT,MTPHI
  47. INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  48. INTEGER IAROTA(NSB)
  49. REAL*8 XPHILB(NSB,NPLSB,NA2,IDIMB)
  50. ENDSEGMENT
  51. SEGMENT MTLIAB
  52. INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
  53. REAL*8 XPALB(NLIAB,NXPALB)
  54. REAL*8 XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  55. ENDSEGMENT
  56. SEGMENT,MTNUM
  57. REAL*8 XDT(NPC1),XTEMPS(NPC1)
  58. ENDSEGMENT
  59. SEGMENT,MPREF
  60. INTEGER IPOREF(NPREF)
  61. ENDSEGMENT
  62. *
  63. *
  64. segment mtbas
  65. integer itbmod,lsstru(np1),nsstru
  66. endsegment
  67.  
  68. c segment local : enregistre les positions d'indices
  69. segment IN2IA(LVAL)
  70. c segment local : calcule l operateur et son inverse
  71. SEGMENT MOP
  72. REAL*8 XOP(NB1,NB1)
  73. INTEGER INVOP(NB1)
  74. REAL*8 XOPM1(NB1,NB1)
  75. ENDSEGMENT
  76.  
  77. LOGICAL L0,L1,RIGIDE,LMODYN,LPLEIN(3)
  78. CHARACTER*4 CMOT,MOINC
  79. CHARACTER*8 TYPRET,CHARRE
  80. CHARACTER*40 MONMOT
  81. REAL*8 GAMMA,BETA
  82. DATA GAMMA /0.5D0/
  83. *
  84. MTKAM = KTKAM
  85. MTPHI = KTPHI
  86. MTLIAB = KTLIAB
  87. MPREF = KPREF
  88. MTNUM = KTNUM
  89. *
  90. NPLB = IBASB(/1)
  91. NSB = INMSB(/1)
  92. NA2 = XPHILB(/3)
  93. IDIMB = XPHILB(/4)
  94. NLIAB = IPALB(/1)
  95. NA1 = XASM(/1)
  96. NB1K = XK(/2)
  97. NB1C = XASM(/2)
  98. NB1M = XM(/2)
  99. NB1 = XOPER(/1)
  100. NOPER= XOPER(/3)
  101. NPREF=IPOREF(/1)
  102. *
  103. IA1 = 0
  104. DEUXPI = 2.D0 * XPI
  105. RIGIDE =.FALSE.
  106. LPLEIN(1)=.FALSE.
  107. LPLEIN(2)=.FALSE.
  108. LPLEIN(3)=.FALSE.
  109. *
  110. ************************************************************************
  111. * Traitement des matrices de variables generalisees K et M
  112. ************************************************************************
  113. *
  114. IF (ITBAS.NE.0 .AND.ITKM.EQ.0.AND.(.NOT.LMODYN)) THEN
  115. CALL ACCTAB(ITBAS,'MOT',IMODE,X0,'SOUSTYPE',L0,IP0,
  116. & 'MOT',I1,X1,MONMOT,L1,IP1)
  117. IF (IERR.NE.0) RETURN
  118. *
  119. *------- Cas ou la base est unique ------------------------------------*
  120. *
  121. IF (MONMOT(1:11).EQ.'BASE_MODALE') THEN
  122. *
  123. * On recupere la base de modes
  124. *
  125. CALL ACCTAB(ITBAS,'MOT',IMODE,X0,'MODES',L0,IP0,
  126. & 'TABLE',I1,X1,' ',L1,IBAS)
  127. IF (IERR.NE.0) RETURN
  128. CALL DYNE26(IBAS,KTKAM,KTLIAB,KTPHI,IA1,1,ICOMP,
  129. & RIGIDE,ITCARA,LMODYN,ITKM)
  130. c IF (RIGIDE) THEN
  131. c RIGIDE =.FALSE.
  132. c DO 80 ILIA =1,NLIAB
  133. c ITYP = IPALB(ILIA,1)
  134. c IF (ITYP.EQ.35.OR.ITYP.EQ.36) THEN
  135. c RIGIDE =.TRUE.
  136. c ENDIF
  137. c 80 CONTINUE
  138. c ENDIF
  139. IF (IERR.NE.0) RETURN
  140. *
  141. *------- Cas ou on a un ensemble de bases -----------------------------*
  142. *
  143. ELSE IF (MONMOT(1:17).EQ.'ENSEMBLE_DE_BASES') THEN
  144. *
  145. * On boucle sur le nombre de bases
  146. *
  147. IT = 0
  148. NPLSB = 0
  149. 10 CONTINUE
  150. TYPRET = ' '
  151. IT = IT + 1
  152. CALL ACCTAB(ITBAS,'ENTIER',IT,X0,' ',L0,IP0,
  153. & TYPRET,I1,X1,CHARRE,L1,ITTBAS)
  154. IF (IERR.NE.0) RETURN
  155. IF (ITTBAS.NE.0) THEN
  156. IF (TYPRET.EQ.'TABLE ') THEN
  157. CALL ACCTAB(ITTBAS,'MOT',IMODE,X0,'MODES',L0,IP0,
  158. & 'TABLE',I1,X1,' ',L1,IBAS)
  159. IF (IERR.NE.0) RETURN
  160. CALL DYNE26(IBAS,KTKAM,KTLIAB,KTPHI,IA1,IT,ICOMP,
  161. & RIGIDE,ITCARA,LMODYN,ITKM)
  162. IF (IERR.NE.0) RETURN
  163. NPLSB = MAX(NPLSB,ICOMP)
  164. GOTO 10
  165. ELSE
  166. CALL ERREUR(491)
  167. RETURN
  168. ENDIF
  169. ENDIF
  170.  
  171. * le segadj n'est plus necessaire
  172. * on le fait dans dyne26
  173. * MP
  174. * SEGADJ,MTPHI
  175. ENDIF
  176. c write(ioimp,*) 'fin du traitement de la base',ITBAS
  177. *
  178. * ELSE IF (ITBAS.NE.0.AND.ITKM.NE.0) THEN
  179. * WRITE(IOIMP,*)
  180. * & 'DYNE : TBAS et TKM coexistent ---> @ implementer.'
  181. * IERR = 2
  182. * RETURN
  183. *
  184. *------- Cas d'une table PASAPAS en entree ----------------------------*
  185. *
  186. ELSE IF (LMODYN) THEN
  187. mtbas = itbas
  188. NPLSB = 0
  189. do isstru = 1,nsstru
  190. CALL DYNE26(itbas,KTKAM,KTLIAB,KTPHI,IA1,isstru,ICOMP,
  191. & RIGIDE,ITCARA,LMODYN,ITKM)
  192. IF (IERR.NE.0) RETURN
  193. NPLSB = MAX(NPLSB,ICOMP)
  194. enddo
  195. *
  196. *------- Cas d'une table 'RAIDEUR_ET_MASSE' en entree -----------------*
  197. *
  198. ELSE IF (ITKM.NE.0) THEN
  199. TYPRET = ' '
  200. CALL ACCTAB(ITKM,'MOT',I0,X0,'RAIDEUR',L0,IP0,
  201. & TYPRET,I1,X1,CHARRE,L1,IRIGI)
  202. IF (IERR.NE.0) RETURN
  203. IF (IRIGI.NE.0 .AND. TYPRET.EQ.'RIGIDITE') THEN
  204. c nature de la matrice : DIAGONALE (defaut) ou PLEINE ...
  205. TYPRET= ' '
  206. CALL ACCTAB(ITKM,'MOT',I0,X0,'NATURE_RAIDEUR',L0,IP0,
  207. & TYPRET,I1,X1,CHARRE,L1,IP1)
  208. IF(TYPRET(1:3).eq.'MOT') THEN
  209. if(iimpi.EQ.333) write(ioimp,*) 'Nature de K =',CHARRE
  210. IF(CHARRE(1:6).eq.'PLEINE') THEN
  211. LPLEIN(1)=.TRUE.
  212. NB1K=NA1
  213. NB1=max(NB1,NB1K)
  214. SEGADJ,MTKAM
  215. ELSEIF(CHARRE(1:6).ne.'DIAGONALE') THEN
  216. write(ioimp,*) 'Nature de K =',CHARRE,' non comprise !'
  217. call erreur(251)
  218. ENDIF
  219. ELSEIF(TYPRET(1:8).NE.' ') THEN
  220. write(ioimp,*) 'NATURE doit etre un MOT (DIAGONALE ou PLEINE)'
  221. write(ioimp,*) 'et pas un ',TYPRET
  222. MOTERR(1:8)='MOT '
  223. call erreur(37)
  224. ENDIF
  225. MRIGID = IRIGI
  226. SEGACT,MRIGID
  227. NRIGI = IRIGEL(/2)
  228. DO 20 I=1,NRIGI
  229. COEF = COERIG(I)
  230. MELEME = IRIGEL(1,I)
  231. DESCR = IRIGEL(3,I)
  232. XMATRI = IRIGEL(4,I)
  233. SEGACT,DESCR,MELEME,XMATRI
  234. NRIG = RE(/3)
  235. LVAL = RE(/1)
  236. DO 30 IRIG=1,NRIG
  237. IF(LPLEIN(1)) segini,IN2IA
  238. c boucle sur les lignes (ddls duals)
  239. DO 35 IN=1,LVAL
  240. INODE=NOELED(IN)
  241. IF(INODE.ne.NOELEP(IN)) THEN
  242. WRITE(IIOMP,*) 'Incoherence entre les inconnues',
  243. & 'primales et duales de la matrice RAIDEUR'
  244. CALL ERREUR(47)
  245. RETURN
  246. ENDIF
  247. NNODE=NUM(INODE,IRIG)
  248. c position de cette inconnue dans IPOREF de MPREF
  249. DO 36 IA=1,NPREF
  250. IF (IPOREF(IA).EQ.NNODE) GOTO 39
  251. 36 CONTINUE
  252. write(ioimp,*) 'D2VTRA: Incoherence entre les ',
  253. & 'points de reference et la matrice RAIDEUR'
  254. CALL ERREUR(504)
  255. 39 CONTINUE
  256. c write(ioimp,*) 'DEVTRA: + noeud dual trouvé en position',IA
  257. IF(LPLEIN(1)) THEN
  258. c on enregistre la position
  259. IN2IA(IN)=IA
  260. c boucle sur les ddl duals JN >= IN (depuis le coin)
  261. DO 37 JN=1,IN
  262. IB = IN2IA(JN)
  263. * Matrice pleine ...
  264. XK(IA,IB) = XK(IA,IB)
  265. & + (RE(IN,JN,IRIG) * COEF)
  266. c attention a ne pas remplir 2 fois la diagonale...
  267. IF(IA.eq.IB) GOTO 37
  268. XK(IB,IA) = XK(IB,IA)
  269. & + (RE(JN,IN,IRIG) * COEF)
  270. 37 CONTINUE
  271. ELSE
  272. * Partie diagonale seulement ...
  273. XK(IA,1) = XK(IA,1) + (RE(IN,IN,IRIG) * COEF)
  274. ENDIF
  275. 35 CONTINUE
  276. IF(LPLEIN(1)) segsup,IN2IA
  277. 30 CONTINUE
  278. SEGDES,XMATRI,MELEME,DESCR
  279. 20 CONTINUE
  280. SEGDES,MRIGID
  281. ELSE
  282. CALL ERREUR(483)
  283. RETURN
  284. ENDIF
  285. *
  286. TYPRET = ' '
  287. CALL ACCTAB(ITKM,'MOT',I0,X0,'MASSE',L0,IP0,
  288. & TYPRET,I1,X1,CHARRE,L1,IMASS)
  289. IF (IERR.NE.0) RETURN
  290. IF (IMASS.NE.0 .AND. TYPRET.EQ.'RIGIDITE') THEN
  291. c nature de la matrice : DIAGONALE (defaut) ou PLEINE ...
  292. TYPRET= ' '
  293. CALL ACCTAB(ITKM,'MOT',I0,X0,'NATURE_MASSE',L0,IP0,
  294. & TYPRET,I1,X1,CHARRE,L1,IP1)
  295. IF(TYPRET(1:3).eq.'MOT') THEN
  296. if(iimpi.EQ.333) write(ioimp,*) 'Nature de M =',CHARRE
  297. IF(CHARRE(1:6).eq.'PLEINE') THEN
  298. LPLEIN(3)=.TRUE.
  299. NB1M=NA1
  300. NB1=max(NB1,NB1M)
  301. NOPER=1
  302. SEGADJ,MTKAM
  303. ELSEIF(CHARRE(1:6).ne.'DIAGONALE') THEN
  304. write(ioimp,*) 'Nature de M =',CHARRE,' non comprise !'
  305. call erreur(251)
  306. ENDIF
  307. ELSEIF(TYPRET(1:8).NE.' ') THEN
  308. write(ioimp,*) 'NATURE doit etre un MOT (DIAGONALE ou PLEINE)'
  309. write(ioimp,*) 'et pas un ',TYPRET
  310. MOTERR(1:8)='MOT '
  311. call erreur(37)
  312. ENDIF
  313. MRIGID = IMASS
  314. SEGACT,MRIGID
  315. NMASS = IRIGEL(/2)
  316. DO 40 I=1,NMASS
  317. COEF = COERIG(I)
  318. MELEME = IRIGEL(1,I)
  319. DESCR = IRIGEL(3,I)
  320. XMATRI = IRIGEL(4,I)
  321. SEGACT,DESCR,MELEME,XMATRI
  322. NRIG = RE(/3)
  323. LVAL = RE(/1)
  324. DO 50 IRIG=1,NRIG
  325. IF(LPLEIN(3)) segini,IN2IA
  326. c boucle sur les lignes (ddls duals)
  327. DO 55 IN=1,LVAL
  328. INODE=NOELED(IN)
  329. IF(INODE.ne.NOELEP(IN)) THEN
  330. WRITE(IIOMP,*) 'Incoherence entre les inconnues',
  331. & 'primales et duales de la matrice MASSE'
  332. CALL ERREUR(47)
  333. RETURN
  334. ENDIF
  335. NNODE=NUM(INODE,IRIG)
  336. c position de cette inconnue dans IPOREF de MPREF
  337. DO 56 IA=1,NPREF
  338. IF (IPOREF(IA).EQ.NNODE) GOTO 59
  339. 56 CONTINUE
  340. write(ioimp,*) 'D2VTRA: Incoherence entre les ',
  341. & 'points de reference et la matrice MASSE'
  342. CALL ERREUR(504)
  343. 59 CONTINUE
  344. IF(LPLEIN(3)) THEN
  345. c on enregistre la position
  346. IN2IA(IN)=IA
  347. c boucle sur les ddl duals JN >= IN (depuis le coin)
  348. DO 57 JN=1,IN
  349. IB = IN2IA(JN)
  350. * Matrice pleine ...
  351. XM(IA,IB) = XM(IA,IB)
  352. & + (RE(IN,JN,IRIG) * COEF)
  353. c attention a ne pas remplir 2 fois la diagonale...
  354. IF(IA.eq.IB) GOTO 57
  355. XM(IB,IA) = XM(IB,IA)
  356. & + (RE(JN,IN,IRIG) * COEF)
  357. 57 CONTINUE
  358. ELSE
  359. * Partie diagonale seulement ...
  360. XM(IA,1) = XM(IA,1) + (RE(IN,IN,IRIG) * COEF)
  361. ENDIF
  362. 55 CONTINUE
  363. IF(LPLEIN(3)) segsup,IN2IA
  364. 50 CONTINUE
  365. SEGDES,XMATRI,MELEME,DESCR
  366. 40 CONTINUE
  367. SEGDES,MRIGID
  368. ELSE
  369. CALL ERREUR(484)
  370. RETURN
  371. ENDIF
  372.  
  373. * traitement de la base modale (necessaire si liaison B)
  374. * --> remplissage de XPHILB
  375. TYPRET = ' '
  376. CALL ACCTAB(ITKM,'MOT',I0,X0,'BASE_MODALE',L0,IP0,
  377. & TYPRET,I1,X1,CHARRE,L1,ITBAS2)
  378. IF(ITBAS2.eq.0) GOTO 599
  379. TYPRET = ' '
  380. CALL ACCTAB(ITBAS2,'MOT',IMODE,X0,'SOUSTYPE',L0,IP0,
  381. & 'MOT',I1,X1,MONMOT,L1,IP1)
  382. IF (IERR.NE.0) RETURN
  383. * Cas ou la base est unique
  384. IF (MONMOT(1:11).EQ.'BASE_MODALE') THEN
  385. * On recupere la base de modes
  386. CALL ACCTAB(ITBAS2,'MOT',IMODE,X0,'MODES',L0,IP0,
  387. & 'TABLE',I1,X1,' ',L1,IBAS2)
  388. IF (IERR.NE.0) RETURN
  389. CALL DYNE26(IBAS2,KTKAM,KTLIAB,KTPHI,IA1,1,ICOMP,
  390. & RIGIDE,ITCARA,LMODYN,ITKM)
  391. * Cas ou la base est un ensemble de bases
  392. ELSEIF(MONMOT(1:17).EQ.'ENSEMBLE_DE_BASES') THEN
  393. IT = 0
  394. NPLSB = 0
  395. 510 CONTINUE
  396. TYPRET = ' '
  397. IT = IT + 1
  398. CALL ACCTAB(ITBAS2,'ENTIER',IT,X0,' ',L0,IP0,
  399. & TYPRET,I1,X1,CHARRE,L1,ITTBAS)
  400. IF (IERR.NE.0) RETURN
  401. IF (ITTBAS.NE.0) THEN
  402. IF (TYPRET.EQ.'TABLE ') THEN
  403. CALL ACCTAB(ITTBAS,'MOT',IMODE,X0,'MODES',L0,IP0,
  404. & 'TABLE',I1,X1,' ',L1,IBAS2)
  405. IF (IERR.NE.0) RETURN
  406. CALL DYNE26(IBAS2,KTKAM,KTLIAB,KTPHI,IA1,IT,ICOMP,
  407. & RIGIDE,ITCARA,LMODYN,ITKM)
  408. IF (IERR.NE.0) RETURN
  409. NPLSB = MAX(NPLSB,ICOMP)
  410. GOTO 510
  411. ELSE
  412. CALL ERREUR(491)
  413. RETURN
  414. ENDIF
  415. ENDIF
  416. ENDIF
  417. 599 CONTINUE
  418. IF (IERR.NE.0) RETURN
  419. ENDIF
  420.  
  421.  
  422. ************************************************************************
  423. * Traitement de la matrice d'amortissement
  424. ************************************************************************
  425. *
  426. IF (ITA.NE.0) THEN
  427. if (lmodyn) then
  428. iamor = ita
  429. typret='RIGIDITE'
  430. else
  431. TYPRET = ' '
  432. CALL ACCTAB(ITA,'MOT',I0,X0,'AMORTISSEMENT',L0,IP0,
  433. & TYPRET,I1,X1,CHARRE,L1,IAMOR)
  434. IF (IERR.NE.0) RETURN
  435. endif
  436. IF (IAMOR.NE.0 .AND. TYPRET.EQ.'RIGIDITE') THEN
  437. c nature de la matrice : DIAGONALE (defaut) ou PLEINE ...
  438. if(.not.lmodyn) then
  439. TYPRET= ' '
  440. CALL ACCTAB(ITA,'MOT',I0,X0,'NATURE',L0,IP0,
  441. & TYPRET,I1,X1,CHARRE,L1,IP1)
  442. IF(TYPRET(1:3).eq.'MOT') THEN
  443. if(iimpi.EQ.333) write(ioimp,*) 'Nature de C =',CHARRE
  444. IF(CHARRE(1:6).eq.'PLEINE') THEN
  445. LPLEIN(2)=.TRUE.
  446. NB1C=NA1
  447. NB1=max(NB1,NB1C)
  448. NOPER=max(1,NOPER)
  449. SEGADJ,MTKAM
  450. ELSEIF(CHARRE(1:6).ne.'DIAGONALE') THEN
  451. write(ioimp,*) 'Nature de C =',CHARRE,' non comprise !'
  452. call erreur(251)
  453. ENDIF
  454. ELSEIF(TYPRET(1:8).NE.' ') THEN
  455. write(ioimp,*) 'NATURE doit etre un MOT (DIAGONALE ou PLEINE)'
  456. write(ioimp,*) 'et pas un ',TYPRET
  457. MOTERR(1:8)='MOT '
  458. call erreur(37)
  459. ENDIF
  460. endif
  461. MRIGID = IAMOR
  462. SEGACT,MRIGID
  463. NAMOR = IRIGEL(/2)
  464. DO 60 I=1,NAMOR
  465. COEF = COERIG(I)
  466. c write(ioimp,*) 'D2VTRA: sous rigidite ',I,'/',NAMOR,COEF
  467. MELEME = IRIGEL(1,I)
  468. DESCR = IRIGEL(3,I)
  469. XMATRI = IRIGEL(4,I)
  470. SEGACT,DESCR,MELEME,XMATRI
  471. NRIG = RE(/3)
  472. LVAL = RE(/1)
  473. DO 70 IRIG=1,NRIG
  474. c write(ioimp,*) 'D2VTRA: + element',IRIG,'/',NRIG
  475. IF(LPLEIN(2)) segini,IN2IA
  476. c boucle sur les lignes (ddls duals)
  477. DO 75 IN=1,LVAL
  478. INODE=NOELED(IN)
  479. IF(INODE.ne.NOELEP(IN)) THEN
  480. WRITE(IIOMP,*) 'Incoherence entre les inconnues',
  481. & 'primales et duales de la matrice AMORTISSEMENT'
  482. CALL ERREUR(47)
  483. RETURN
  484. ENDIF
  485. NNODE=NUM(INODE,IRIG)
  486. c write(ioimp,*) 'D2VTRA: + noeud dual',IN,'/',LVAL,' #',NNODE
  487. c position de cette inconnue dans IPOREF de MPREF
  488. DO 76 IA=1,NPREF
  489. IF (IPOREF(IA).EQ.NNODE) GOTO 79
  490. 76 CONTINUE
  491. write(ioimp,*) 'D2VTRA: Incoherence entre les ',
  492. & 'points de reference et la matrice AMORTISSEMENT'
  493. CALL ERREUR(504)
  494. 79 CONTINUE
  495. c write(ioimp,*) 'D2VTRA: + noeud dual trouvé en position',IA
  496. IF(LPLEIN(2)) THEN
  497. c on enregistre la position
  498. IN2IA(IN)=IA
  499. c boucle sur les ddl duals JN >= IN (depuis le coin)
  500. DO 77 JN=1,IN
  501. IB = IN2IA(JN)
  502. * Matrice pleine ...
  503. XASM(IA,IB) = XASM(IA,IB)
  504. & + (RE(IN,JN,IRIG) * COEF)
  505. c attention a ne pas remplir 2 fois la diagonale...
  506. IF(IA.eq.IB) GOTO 77
  507. XASM(IB,IA) = XASM(IB,IA)
  508. & + (RE(JN,IN,IRIG) * COEF)
  509. 77 CONTINUE
  510. ELSE
  511. * Partie diagonale seulement ...
  512. XASM(IA,1) = XASM(IA,1) + (RE(IN,IN,IRIG) * COEF)
  513. ENDIF
  514. 75 CONTINUE
  515. IF(LPLEIN(2)) segsup,IN2IA
  516. 70 CONTINUE
  517. SEGDES,XMATRI,MELEME,DESCR
  518. 60 CONTINUE
  519. SEGDES,MRIGID
  520. ELSE
  521. CALL ERREUR(485)
  522. RETURN
  523. ENDIF
  524. ENDIF
  525.  
  526. ************************************************************************
  527. * on calcule 1 operateur = inverse de :
  528. * IALGO=2 : XOP = M + 0.5dt*C
  529. * IALGO=3 : XOP = M + 0.5dt*C
  530. * IALGO=4 : XOP = M + 0.5dt*C
  531. ************************************************************************
  532.  
  533. c IF(LPLEIN(2).OR.LPLEIN(3)) THEN
  534. IF(NOPER.GT.0) THEN
  535. SEGINI,MOP
  536. c le pas de temps est constant (seule possibilite aujourd'hui)
  537. B3 = GAMMA*xdt(1)
  538. c ...C pleine, M pleine
  539. IF(LPLEIN(2).AND.LPLEIN(3)) THEN
  540. DO 61 IA=1,NA1
  541. DO 61 IB=1,NA1
  542. XOP(IA,IB) = XM(IA,IB) + B3*XASM(IA,IB)
  543. 61 CONTINUE
  544. c ...C pleine, M diago
  545. ELSEIF(LPLEIN(2)) THEN
  546. DO 62 IA=1,NA1
  547. XOP(IA,IA) = XM(IA,1)
  548. DO 62 IB=1,NA1
  549. XOP(IA,IB) = XOP(IA,IB) + B3*XASM(IA,IB)
  550. 62 CONTINUE
  551. c ...C diago, M pleine
  552. ELSEIF(LPLEIN(3)) THEN
  553. DO 63 IA=1,NA1
  554. XOP(IA,IA) = B3*XASM(IA,1)
  555. DO 63 IB=1,NA1
  556. XOP(IA,IB) = XOP(IA,IB) + XM(IA,IB)
  557. 63 CONTINUE
  558. ELSEIF(.NOT.LPLEIN(1).OR.IALGO.LE.2) THEN
  559. WRITE(IOIMP,*) 'D2VTRA: option PLEINE incoherente !'
  560. CALL ERREUR(5)
  561. RETURN
  562. ENDIF
  563. c Ajout de K si implicite
  564. IF(IALGO.GE.3) THEN
  565. B4 = BETA*(xdt(1)**2)
  566. IF(LPLEIN(1)) THEN
  567. DO 64 IA=1,NA1
  568. DO 64 IB=1,NA1
  569. XOP(IA,IB) = XOP(IA,IB) + B4*XK(IA,IB)
  570. 64 CONTINUE
  571. ELSE
  572. DO 65 IA=1,NA1
  573. XOP(IA,IA) = XOP(IA,IA) + B4*XK(IA,1)
  574. 65 CONTINUE
  575. ENDIF
  576. ENDIF
  577. * -Inversion de l'operateur
  578. CALL IVMAT(NA1,XOP,INVOP,XOPM1,DETOP,0,IRET)
  579. IF(IRET.ne.0) RETURN
  580. * rem : ici IVMAT, mais il existe aussi INVER, INVER1 ...
  581. DO 69 IA=1,NA1
  582. DO 69 IB=1,NB1
  583. XOPER(IA,IB,1)=XOPM1(IA,IB)
  584. 69 CONTINUE
  585. SEGSUP,MOP
  586. ENDIF
  587. *
  588. IF (IIMPI.EQ.333) THEN
  589. WRITE(IOIMP,*)' segment MTPHI'
  590. WRITE(IOIMP,*)'D2VTRA : valeur de NPLB :',IBASB(/1)
  591. WRITE(IOIMP,*)'D2VTRA : valeur de NSB :',XPHILB(/1)
  592. WRITE(IOIMP,*)'D2VTRA : valeur de NPLSB :',XPHILB(/2)
  593. WRITE(IOIMP,*)'D2VTRA : valeur de NA2 :',XPHILB(/3)
  594. WRITE(IOIMP,*)'D2VTRA : valeur de IDIMB :',XPHILB(/4)
  595. if(NB1K.gt.1) then
  596. do iou=1,NA1
  597. WRITE(IOIMP,*) 'XK=',(XK(iou,jou),jou=1,NB1K)
  598. enddo
  599. else
  600. do iou=1,NA1
  601. WRITE(IOIMP,*) 'XK(',iou,',1)=',XK(iou,1)
  602. enddo
  603. endif
  604. if(NB1C.gt.1) then
  605. do iou=1,NA1
  606. WRITE(IOIMP,*) 'XASM=',(XASM(iou,jou),jou=1,NB1C)
  607. enddo
  608. else
  609. do iou=1,NA1
  610. WRITE(IOIMP,*) 'XASM(',iou,',1)=',XASM(iou,1)
  611. enddo
  612. endif
  613. if(NB1M.gt.1) then
  614. do iou=1,NA1
  615. WRITE(IOIMP,*) 'XM=',(XM(iou,jou),jou=1,NB1M)
  616. enddo
  617. else
  618. do iou=1,NA1
  619. WRITE(IOIMP,*) 'XM(',iou,',1)=',XM(iou,1)
  620. enddo
  621. endif
  622. if(NOPER.ge.1) then
  623. do iop=1,NOPER
  624. do iou=1,NB1
  625. WRITE(IOIMP,*) 'XOPER(',iou,',:,',iop,')=',
  626. & (XOPER(iou,jou,iop),jou=1,NB1)
  627. enddo
  628. enddo
  629. endif
  630. ENDIF
  631. *
  632. END
  633.  
  634.  
  635.  
  636.  
  637.  
  638.  
  639.  
  640.  
  641.  
  642.  
  643.  
  644.  
  645.  
  646.  
  647.  
  648.  
  649.  

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