Télécharger devtra.eso

Retour à la liste

Numérotation des lignes :

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

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