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

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