Télécharger devini.eso

Retour à la liste

Numérotation des lignes :

  1. C DEVINI SOURCE BP208322 15/07/22 21:15:20 8586
  2. SUBROUTINE DEVINI(ITINIT,KTKAM,KTQ,KTFEX,KTPAS,KTNUM,KTLIAA,
  3. & KTLIAB,KTPHI,KCPR,KOCLFA,KOCLB1,REPRIS,
  4. & RIGIDE,lmodyn)
  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. * Initialisation de l'algorithme ou reprise de calcul. *
  13. * *
  14. * Parametres: *
  15. * *
  16. * e ITINIT Table contenant les conditions initiales ou les *
  17. * champs necessaires a la reprise du calcul *
  18. * es KTKAM Segment contenant les vecteurs XK, XASM et XM *
  19. * es KTQ Segment des variables de mouvement generalisees *
  20. * (et des travaux)
  21. * es KTFEX Segment contenant les chargements libres *
  22. * es KTPAS Segment des variables au cours d'un pas de temps *
  23. * es KTNUM Segment contenant les parametres numeriques *
  24. * e KTLIAA Segment descriptif des liaisons en base A *
  25. * e KTLIAB Segment descriptif des liaisons en base B *
  26. * e KTPHI Segment contenant les deformees modales *
  27. * e KCPR Segment des points *
  28. * - KOCLFA Segment contenant les tableaux locaux de la subroutine *
  29. * DEVLFA *
  30. * - KOCLB1 Segment contenant les tableaux locaux de la subroutine *
  31. * DEVLB1 *
  32. * e REPRIS Vrai si reprise de calcul *
  33. * e RIGIDE Vrai si corps rigide
  34. * *
  35. * Auteur, date de creation: *
  36. * *
  37. * Denis ROBERT-MOUGIN, le 25 mai 1989. *
  38. * *
  39. *--------------------------------------------------------------------*
  40. -INC CCOPTIO
  41. -INC CCNOYAU
  42. -INC SMCOORD
  43. -INC SMTABLE
  44. -INC CCASSIS
  45. *
  46. SEGMENT,MTQ
  47. REAL*8 Q1(NA1,4),Q2(NA1,4),Q3(NA1,4)
  48. REAL*8 WEXT(NA1,2),WINT(NA1,2)
  49. ENDSEGMENT
  50. SEGMENT,MTFEX
  51. REAL*8 FEXA(NPFEXA,NPC1,2)
  52. REAL*8 FEXPSM(NPLB,NPC1,2,IDIMB)
  53. REAL*8 FTEXB(NPLB,NPC1,2,IDIM)
  54. INTEGER IFEXA(NPFEXA),IFEXB(NPFEXB)
  55. ENDSEGMENT
  56. SEGMENT,MTPAS
  57. REAL*8 FTOTA(NA1,4),FTOTB(NPLB,IDIMB),FTOTBA(NA1)
  58. REAL*8 XPTB(NPLB,4,IDIMB),FINERT(NA1,4)
  59. REAL*8 XVALA(NLIAA,4,NTVAR),XVALB(NLIAB,4,NTVAR)
  60. REAL*8 FEXB(NPLB,2,IDIM),XCHPFB(2,NLIAB,4,NPLB)
  61. ENDSEGMENT
  62. SEGMENT,MTNUM
  63. REAL*8 XDT(NPC1),XTEMPS(NPC1)
  64. ENDSEGMENT
  65. SEGMENT,MTKAM
  66. REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
  67. REAL*8 XOPER(NB1,NB1,NOPER)
  68. ENDSEGMENT
  69. SEGMENT,MTLIAA
  70. INTEGER IPALA(NLIAA,NIPALA),IPLIA(NLIAA,NPLAA),JPLIA(NPLA)
  71. REAL*8 XPALA(NLIAA,NXPALA)
  72. ENDSEGMENT
  73. SEGMENT,MTLIAB
  74. INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
  75. REAL*8 XPALB(NLIAB,NXPALB)
  76. REAL*8 XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  77. ENDSEGMENT
  78. SEGMENT,MTPHI
  79. INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  80. INTEGER IAROTA(NSB)
  81. REAL*8 XPHILB(NSB,NPLSB,NA2,IDIMB)
  82. ENDSEGMENT
  83. SEGMENT,ICPR(XCOOR(/1)/(IDIM+1))
  84.  
  85. * Segment "local" pour DEVLFA ...
  86. SEGMENT,LOCLFA
  87. REAL*8 FTEST(NA1,4),FTOTA0(NA1,4)
  88. ENDSEGMENT
  89. * Segment "local" pour DEVLB1 ...
  90. SEGMENT,LOCLB1
  91. REAL*8 FTEST2(NPLB,6),FTOTB0(NPLB,6)
  92. ENDSEGMENT
  93. segment mwinit
  94. integer jpdep,jpvit,jrepr
  95. endsegment
  96. c segment local
  97. SEGMENT MAMOR
  98. REAL*8 FAMOR(NA1,4)
  99. REAL*8 FTEMP(NA1),QTEMP(NA1)
  100. ENDSEGMENT
  101. c segment local : calcule un operateur et son inverse
  102. SEGMENT MOP
  103. REAL*8 XOP(NB1,NB1)
  104. INTEGER INVOP(NB1)
  105. REAL*8 XOPM1(NB1,NB1)
  106. ENDSEGMENT
  107.  
  108. *
  109. LOGICAL L0,L1,REPRIS,RIGIDE,lmodyn
  110. CHARACTER*8 TYPRET,TYPIND,CHARRE
  111. CHARACTER*(19) CHAI1
  112. ILC1 = 19
  113. DATA CHAI1 /'VARIABLES_LIAISON_A'/
  114. *
  115. MTFEX = KTFEX
  116. MTPAS = KTPAS
  117. MTNUM = KTNUM
  118. MTKAM = KTKAM
  119. MTQ = KTQ
  120. MTLIAA = KTLIAA
  121. MTLIAB = KTLIAB
  122. MTPHI = KTPHI
  123. LOCLFA = KOCLFA
  124. LOCLB1 = KOCLB1
  125. IDEPL = 0
  126. IVITE = 0
  127. ITABL = 0
  128. ICH1 = 0
  129. ICH2 = 0
  130. ICH3 = 0
  131. ICH4 = 0
  132. ICH5 = 0
  133. ICH6 = 0
  134. NA1 = Q1(/1)
  135. NB1K = XK(/2)
  136. NB1C = XASM(/2)
  137. NB1M = XM(/2)
  138. NB1 = XOPER(/1)
  139. NLIAA = IPALA(/1)
  140. NLIAB = IPALB(/1)
  141. NPLB = JPLIB(/1)
  142. NSB = XPHILB(/1)
  143. NPLSB = XPHILB(/2)
  144. NA2 = XPHILB(/3)
  145. IDIMB = XPHILB(/4)
  146. NPFEXA = FEXA(/1)
  147. NPC1 = FEXPSM(/2)
  148. NIP = XABSCI(/2)
  149. IERRD = 0
  150. *
  151. * S'agit-il d'une initialisation ou d'une reprise de calcul ?
  152. *
  153. IF ( REPRIS ) THEN
  154. IF (IIMPI.EQ.333) THEN
  155. WRITE(IOIMP,*)'DEVINI : ---> reprise de calcul'
  156. ENDIF
  157. if (lmodyn) then
  158. mwinit = itinit
  159. segact mwinit
  160. ITABL = jrepr
  161. else
  162. CALL ACCTAB(ITINIT,'MOT',I0,X0,'REPRISE',L0,IP0,
  163. & 'TABLE',I1,X1,' ',L1,ITABL)
  164. endif
  165. * . .
  166. * Reprise du calcul, on remplit: Q Q Q Q
  167. * i,-1/2 i,-1/2 i,0 i,0
  168. *
  169. CALL ACCTAB(ITABL,'MOT',I0,X0,'DEPLACEMENT',L0,IP0,
  170. & 'CHPOINT',I1,X1,' ',L1,ICH1)
  171. IF (IERR.NE.0) RETURN
  172. *
  173. CALL ACCTAB(ITABL,'MOT',I0,X0,'VITESSE',L0,IP0,
  174. & 'CHPOINT',I1,X1,' ',L1,ICH2)
  175. IF (IERR.NE.0) RETURN
  176. *
  177. CALL ACCTAB(ITABL,'MOT',I0,X0,'DEPLACEMENT_1/2',L0,IP0,
  178. & 'CHPOINT',I1,X1,' ',L1,ICH3)
  179. IF (IERR.NE.0) RETURN
  180. *
  181. CALL ACCTAB(ITABL,'MOT',I0,X0,'VITESSE_1/2',L0,IP0,
  182. & 'CHPOINT',I1,X1,' ',L1,ICH4)
  183. IF (IERR.NE.0) RETURN
  184. *
  185. CALL ACCTAB(ITABL,'MOT',I0,X0,'FORCES_BASE_A',L0,IP0,
  186. & 'CHPOINT',I1,X1,' ',L1,ICH5)
  187. IF (IERR.NE.0) RETURN
  188. *
  189. CALL ACCTAB(ITABL,'MOT',I0,X0,'FORCES_BASE_A_1/2',L0,IP0,
  190. & 'CHPOINT',I1,X1,' ',L1,ICH6)
  191. IF (IERR.NE.0) RETURN
  192. *
  193. CALL ACCTAB(ITABL,'MOT',I0,X0,'ACCELERATION',L0,IP0,
  194. & 'CHPOINT',I1,X1,' ',L1,ICH7)
  195. IF (IERR.NE.0) RETURN
  196. *
  197. CALL ACCTAB(ITABL,'MOT',I0,X0,'ACCELERATION_1/2',L0,IP0,
  198. & 'CHPOINT',I1,X1,' ',L1,ICH8)
  199. IF (IERR.NE.0) RETURN
  200. *
  201. CALL ACCTAB(ITABL,'MOT',I0,X0,'TRAVAIL_EXTERIEUR',L0,IP0,
  202. & 'CHPOINT',I1,X1,' ',L1,ICH9)
  203. IF (IERR.NE.0) RETURN
  204. *
  205. CALL ACCTAB(ITABL,'MOT',I0,X0,'TRAVAIL_INTERIEUR',L0,IP0,
  206. & 'CHPOINT',I1,X1,' ',L1,ICH10)
  207. IF (IERR.NE.0) RETURN
  208. *
  209. IF (ICH1.NE.0) THEN
  210. CALL DYNE18(ICH1,KTQ,1,3,KCPR)
  211. ELSE
  212. CALL ERREUR(487)
  213. RETURN
  214. ENDIF
  215. *
  216. IF (ICH2.NE.0) THEN
  217. CALL DYNE18(ICH2,KTQ,2,3,KCPR)
  218. ELSE
  219. CALL ERREUR(487)
  220. RETURN
  221. ENDIF
  222. *
  223. IF (ICH3.NE.0) THEN
  224. CALL DYNE18(ICH3,KTQ,1,4,KCPR)
  225. ELSE
  226. CALL ERREUR(487)
  227. RETURN
  228. ENDIF
  229. *
  230. IF (ICH4.NE.0) THEN
  231. CALL DYNE18(ICH4,KTQ,2,4,KCPR)
  232. ELSE
  233. CALL ERREUR(487)
  234. RETURN
  235. ENDIF
  236. *
  237. IF (ICH5.NE.0) THEN
  238. CALL DYNE23(ICH5,KTPAS,3)
  239. ELSE
  240. CALL ERREUR(487)
  241. RETURN
  242. ENDIF
  243. *
  244. IF (ICH6.NE.0) THEN
  245. CALL DYNE23(ICH6,KTPAS,4)
  246. ELSE
  247. CALL ERREUR(487)
  248. RETURN
  249. ENDIF
  250. *
  251. IF (ICH7.NE.0) THEN
  252. CALL DYNE18(ICH7,KTQ,3,3,KCPR)
  253. ELSE
  254. CALL ERREUR(487)
  255. RETURN
  256. ENDIF
  257. *
  258. IF (ICH8.NE.0) THEN
  259. CALL DYNE18(ICH8,KTQ,3,4,KCPR)
  260. ELSE
  261. CALL ERREUR(487)
  262. RETURN
  263. ENDIF
  264. *
  265. IF (ICH9.NE.0) THEN
  266. CALL DYNE18(ICH9,KTQ,4,1,KCPR)
  267. ENDIF
  268. IF (ICH10.NE.0) THEN
  269. CALL DYNE18(ICH10,KTQ,5,1,KCPR)
  270. ENDIF
  271. *
  272. *
  273. IF (NLIAA.NE.0) THEN
  274. *
  275. * l'indice VARIABLES_LIAISON_A n'existe que pour
  276. * les liaisons POLYNOMIALES
  277. *
  278. IPOLY = 0
  279. MTABLE = ITABL
  280. SEGACT,MTABLE
  281. NIND1 = MLOTAB
  282. if(nbesc.ne.0)segact ipiloc
  283. DO 100 I=1,NIND1
  284. TYPIND = MTABTI(I)
  285. IF (TYPIND.EQ.'MOT ') THEN
  286. IP = MTABII(I)
  287. ID = IPCHAR(IP)
  288. IFI = IPCHAR(IP+1)
  289. IL1 = IFI - ID
  290. IF (IL1.EQ.ILC1) THEN
  291. IF (CHAI1.EQ.ICHARA(ID:IFI-1)) THEN
  292. IPOLY = 1
  293. ENDIF
  294. ENDIF
  295. ENDIF
  296. 100 CONTINUE
  297. if(nbesc.ne.0)segdes ipiloc
  298. SEGDES,MTABLE
  299. IF (IIMPI.EQ.333) THEN
  300. WRITE(IOIMP,*) 'DEVINI : IPOLY = ',IPOLY
  301. ENDIF
  302. IF (IPOLY.NE.0) THEN
  303. CALL ACCTAB(ITABL,'MOT',I0,X0,'VARIABLES_LIAISON_A',
  304. & L0,IP0,'TABLE',I1,X1,' ',L1,ITREFR)
  305. CALL DYNA14(ITREFR,KTLIAA)
  306. IF (IERR.NE.0) RETURN
  307. ENDIF
  308. ENDIF
  309. IF (NLIAB.NE.0) THEN
  310. CALL DEVRCO(Q1,NA1,XPTB,NPLB,XPHILB,NSB,NPLSB,NA2,IDIMB,
  311. & IBASB,IPLSB,INMSB,IORSB,3,IAROTA)
  312. CALL ACCTAB(ITABL,'MOT',I0,X0,'VARIABLES_LIAISON_B',
  313. & L0,IP0,'TABLE',I1,X1,' ',L1,ITREFR)
  314. IF (IERR.NE.0) RETURN
  315. CALL DYNE14(ITREFR,KTLIAB)
  316. IF (IERR.NE.0) RETURN
  317. ENDIF
  318. *
  319. ELSE
  320. *
  321. * Chpoints initiaux de deplacement et de vitesse:
  322. *
  323. IF (ITINIT.NE.0) THEN
  324. if (lmodyn) then
  325. mwinit = itinit
  326. segact mwinit
  327. idepl = jpdep
  328. ivite = jpvit
  329. else
  330. TYPRET=' '
  331. CALL ACCTAB(ITINIT,'MOT',I0,X0,'DEPLACEMENT',L0,IP0,
  332. & TYPRET,I1,X1,CHARRE,L1,IDEPL)
  333. *
  334. TYPRET=' '
  335. CALL ACCTAB(ITINIT,'MOT',I0,X0,'VITESSE',L0,IP0,
  336. & TYPRET,I1,X1,CHARRE,L1,IVITE)
  337. endif
  338. ENDIF
  339. *
  340. * Mise des valeurs initiales au pas m (indice 3)
  341. *
  342. IF (IDEPL.NE.0) THEN
  343. CALL DYNE18(IDEPL,KTQ,1,3,KCPR)
  344. MTQ = KTQ
  345. ENDIF
  346. IVINIT = 0
  347. IF (IVITE.NE.0) THEN
  348. CALL DYNE18(IVITE,KTQ,2,3,KCPR)
  349. IVINIT = 1
  350. ENDIF
  351. *
  352. * Ajout des forces de liaison
  353. *
  354. IF (NLIAA.NE.0) THEN
  355. CALL DEVLFA(Q1,Q2,FTOTA,NA1,IPALA,IPLIA,XPALA,XVALA,
  356. & NLIAA,XDT,1,3,FINERT,IVINIT,FTEST,FTOTA0)
  357. ENDIF
  358. IF (NLIAB.NE.0) THEN
  359. IF (RIGIDE) THEN
  360. DO 13 IP=1,NPLB
  361. DO 15 ID=1,IDIM
  362. FEXB(IP,1,ID) = FTEXB(IP,1,1,ID)
  363. 15 CONTINUE
  364. 13 CONTINUE
  365. ENDIF
  366. c SEGMENT,MTPHI
  367. c INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  368. c INTEGER IAROTA(NSB)
  369. c REAL*8 XPHILB(NSB,NPLSB,NA2,IDIMB)
  370. c ENDSEGMENT
  371. CALL DEVLF2(Q1,Q2,FTOTA,NA1,IPALB,IPLIB,XPALB,XVALB,
  372. & NLIAB,XPHILB,JPLIB,NPLB,IDIMB,FTOTB,FTOTBA,XPTB,
  373. & XDT,1,IBASB,IPLSB,INMSB,IORSB,NSB,NPLSB,NA2,3,
  374. & FEXPSM,NPC1,IERRD,FTEST2,FTOTB0,XABSCI,XORDON,
  375. & NIP,IAROTA,RIGIDE,FEXB,XCHPFB)
  376. IF (IERRD.NE.0) THEN
  377. IF (IERRD.EQ.1) CALL ERREUR(528)
  378. RETURN
  379. ENDIF
  380. ENDIF
  381. IF (NLIAA.NE.0 .OR. NLIAB.NE.0) THEN
  382. DO 5 I=1,NA1
  383. FINERT(I,4) = FINERT(I,3)
  384. 5 CONTINUE
  385. * end do
  386. ENDIF
  387. *
  388. * Calcul des deplacements et des vitesses au pas -1/2 (indice 4)
  389. *
  390. PDT = XDT(1)
  391. PDTS2 = PDT / 2.D0
  392. AUX2 = PDT * PDT / 8.D0
  393. IF (IIMPI.EQ.333) THEN
  394. DO 91 I=1,NA1
  395. WRITE(IOIMP,*)'DEVINI :'
  396. WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',3) =',FTOTA(I,3)
  397. WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',4) =',FTOTA(I,4)
  398. WRITE(IOIMP,*)'DEVINI : Fexa (',I,',1) =',Fexa(I,1)
  399. WRITE(IOIMP,*)'DEVINI : Fexa (',I,',2) =',Fexa(I,2)
  400. 91 CONTINUE
  401. * end do
  402. ENDIF
  403.  
  404. c Cas selon nature (pleine ou diag) des matrices K,C et M :
  405.  
  406. c------- Cas K,C ou M pleine -----------------------------------
  407. IF(NB1K.GT.1.OR.NB1C.GT.1.OR.NB1M.GT.1) THEN
  408. IF (IIMPI.EQ.333) WRITE(IOIMP,*)'pleine K,C,M',NB1K,NB1C,NB1M
  409.  
  410. c F_0 = F^liaison_0 + F^ext_0 - K * Q_0
  411. DO 61 I = 1,NA1
  412. FTOTA(I,3) = FTOTA(I,3) + FEXA(I,1,1)
  413. 61 CONTINUE
  414. CALL DEVLK0(Q1,XK,FTOTA,NA1,NB1K,3)
  415.  
  416. c F^AMOR_0 = C * \dot Q_0
  417. SEGINI,MAMOR
  418. CALL DEVLC0(Q2,XASM,FAMOR,NA1,NB1C,3)
  419.  
  420. c \ddot Q_0 = M-1 * (F_0 - F^AMOR_0)
  421. CALL DEVLM0(Q3,XM,XOPER,FTOTA,FAMOR,NA1,NB1,3,NB1M,3)
  422.  
  423. c Q_-1/2 = Q_0 - dt/2 * \dot Q_0 + dt^2/8 * \ddot Q_0
  424. DO 62 I = 1,NA1
  425. Q1(I,4) = Q1(I,3) - PDTS2*Q2(I,3) + AUX2*Q3(I,3)
  426. 62 CONTINUE
  427.  
  428. c F_-1/2 = F^liaison_-1/2 + F^ext_-1/2 - K * Q_-1/2
  429. DO 63 I = 1,NA1
  430. FTOTA(I,4) = FTOTA(I,4) + FEXA(I,1,2)
  431. 63 CONTINUE
  432. CALL DEVLK0(Q1,XK,FTOTA,NA1,NB1K,4)
  433.  
  434. c \dot Q_-1/2 = [4I+dt*M-1*C]-1
  435. c * ( 4 \dot Q_0 - dt*M-1*(F_-1/2+F_0-F^AMOR_0) )
  436. DO 64 I = 1,NA1
  437. FTEMP(I) = FTOTA(I,4) + FTOTA(I,3) - FAMOR(I,3)
  438. 64 CONTINUE
  439. CALL DEVLM1(QTEMP,XM,XOPER,FTEMP,NA1,NB1,3,NB1M)
  440. DO 65 I = 1,NA1
  441. FTEMP(I) = 4.D0*Q2(I,3) - PDT*QTEMP(I)
  442. 65 CONTINUE
  443. IF(NB1C.GT.1.OR.NB1M.GT.1) THEN
  444. * initialisation de l'operateur (local) et de son inverse
  445. * XOP = 4I - dt*M^-1*C
  446. SEGINI,MOP
  447. c ...C pleine, M pleine
  448. IF(NB1C.GT.1.AND.NB1M.GT.1) THEN
  449. DO 662 IA=1,NA1
  450. XOP(IA,IA) = 4.D0
  451. DO 662 IB=1,NA1
  452. AUX = 0.D0
  453. DO 663 IC=1,NA1
  454. AUX = AUX + XOPER(IA,IC,3)*XASM(IC,IB)
  455. 663 CONTINUE
  456. XOP(IA,IB) = XOP(IA,IB) - PDT*AUX
  457. 662 CONTINUE
  458. c ...C pleine, M diago
  459. ELSEIF(NB1C.GT.1) THEN
  460. DO 664 IA=1,NA1
  461. XOP(IA,IA) = 4.D0
  462. DO 664 IB=1,NA1
  463. XOP(IA,IB) = XOP(IA,IB) - PDT*XASM(IA,IB)/XM(IA,1)
  464. 664 CONTINUE
  465. c ...C diago, M pleine
  466. ELSEIF(NB1M.GT.1) THEN
  467. DO 665 IA=1,NA1
  468. XOP(IA,IA) = 4.D0
  469. DO 665 IB=1,NA1
  470. XOP(IA,IB) = XOP(IA,IB) - PDT*XOPM1(IA,IB)*XASM(IB,1)
  471. 665 CONTINUE
  472. ENDIF
  473. CALL IVMAT(NA1,XOP,INVOP,XOPM1,DETOP,0,IRET)
  474. IF (IIMPI.EQ.333) THEN
  475. WRITE(IOIMP,*) 'FTEMP(:)=',(FTEMP(jou),jou=1,NA1)
  476. do iou=1,NB1
  477. WRITE(IOIMP,*) 'XOPM1(',iou,',:)=',
  478. & (XOPM1(iou,jou),jou=1,NB1)
  479. enddo
  480. ENDIF
  481. c \dot Q_-1/2 = [4I-dt*M-1*C]-1 * FTEMP
  482. DO 67 IA = 1,NA1
  483. Q2(IA,4) = 0.D0
  484. DO 67 IB = 1,NA1
  485. Q2(IA,4) = Q2(IA,4) + XOPM1(IA,IB)*FTEMP(IB)
  486. 67 CONTINUE
  487. SEGSUP,MOP
  488. ELSE
  489. DO 68 IA = 1,NA1
  490. Q2(IA,4) = FTEMP(IA) / (4.D0 - (PDT*XASM(IA,1)/XM(IA,1)))
  491. 68 CONTINUE
  492. ENDIF
  493.  
  494. c \ddot Q_-1/2 = M-1 * ( F_-1/2 - F^AMOR_-1/2)
  495. CALL DEVLC0(Q2,XASM,FAMOR,NA1,NB1C,4)
  496. CALL DEVLM0(Q3,XM,XOPER,FTOTA,FAMOR,NA1,NB1,3,NB1M,4)
  497.  
  498. SEGSUP,MAMOR
  499.  
  500. c------- Cas K,C et M diagonal -----------------------------------
  501. ELSE
  502.  
  503. DO 10 I = 1,NA1
  504. UNSM3 = 1.D0 / ( XM(I,1) - FINERT(I,3) )
  505. UNSM4 = 1.D0 / ( XM(I,1) - FINERT(I,4) )
  506. *
  507. * G = F - K Q
  508. * i,0 i,0 i i,0
  509. *
  510. FTOTA(I,3) = FTOTA(I,3) + FEXA(I,1,1) - (XK(I,1)*Q1(I,3))
  511. AUX1 = (FTOTA(I,3)*UNSM3) - (XASM(I,1)*UNSM3*Q2(I,3))
  512. AUX3 = Q1(I,3) - ( PDTS2 * Q2(I,3) )
  513. *
  514. * . 2 .
  515. * Q = Q - h/2 Q + h /8 ( G - A Q )
  516. * i,-1/2 i,0 i,0 i,0 i i,0
  517. *
  518. Q1(I,4) = AUX3 + ( AUX2 * AUX1 )
  519. *
  520. * G = F - K Q
  521. * i,-1/2 i,-1/2 i i,-1/2
  522. *
  523. FTOTA(I,4) = FTOTA(I,4) + FEXA(I,1,2) - (XK(I,1)*Q1(I,4))
  524. AUX4 = 4.D0 + ( XASM(I,1) * PDT * UNSM3 )
  525. AUX5 = ( FTOTA(I,4) * UNSM4 ) + ( FTOTA(I,3) * UNSM3 )
  526. AUX6 = 1.D0 / ( 4.D0 - ( XASM(I,1) * PDT * UNSM3 ) )
  527. *
  528. * . 1 .
  529. * Q = ________ ( ( 4 + A h ) Q - h ( G + G ) )
  530. * i,-1/2 4 - A h i i,0 i,-1/2 i,0
  531. * i
  532. *
  533. Q2(I,4) = AUX6 * ( (AUX4*Q2(I,3)) - (PDT*AUX5) )
  534. *
  535. * " .
  536. * Q = ( G - A Q ) / M
  537. * i i i i i
  538. *
  539. Q3(I,3) = (FTOTA(I,3)*UNSM3) - (XASM(I,1)*UNSM3*Q2(I,3))
  540. Q3(I,4) = (FTOTA(I,4)*UNSM4) - (XASM(I,1)*UNSM4*Q2(I,4))
  541. 10 CONTINUE
  542.  
  543. ENDIF
  544. *
  545. ENDIF
  546. IF (IIMPI.EQ.333) THEN
  547. DO 30 I=1,NA1
  548. WRITE(IOIMP,*)'DEVINI :'
  549. WRITE(IOIMP,*)'DEVINI : Q1(',I,',3) =',Q1(I,3)
  550. WRITE(IOIMP,*)'DEVINI : Q2(',I,',3) =',Q2(I,3)
  551. WRITE(IOIMP,*)'DEVINI : Q3(',I,',3) =',Q3(I,3)
  552. WRITE(IOIMP,*)'DEVINI :'
  553. WRITE(IOIMP,*)'DEVINI : Q1(',I,',4) =',Q1(I,4)
  554. WRITE(IOIMP,*)'DEVINI : Q2(',I,',4) =',Q2(I,4)
  555. WRITE(IOIMP,*)'DEVINI : Q3(',I,',4) =',Q3(I,4)
  556. 30 CONTINUE
  557. * end do
  558. ENDIF
  559. IF (IIMPI.EQ.333) THEN
  560. DO 31 I=1,NA1
  561. WRITE(IOIMP,*)'DEVINI :'
  562. WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',3) =',FTOTA(I,3)
  563. WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',4) =',FTOTA(I,4)
  564. 31 CONTINUE
  565. * end do
  566. ENDIF
  567.  
  568. *
  569. ICPR = KCPR
  570. SEGSUP,ICPR
  571. *
  572. END
  573.  
  574.  
  575.  
  576.  
  577.  
  578.  
  579.  
  580.  

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