Télécharger devini.eso

Retour à la liste

Numérotation des lignes :

  1. C DEVINI SOURCE BP208322 17/07/10 21:15:04 9488
  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. * et liaison COUPLAGE_DEPLACEMENT avec CONVOLUTION
  278. *
  279. IPOLY = 0
  280. MTABLE = ITABL
  281. SEGACT,MTABLE
  282. NIND1 = MLOTAB
  283. if(nbesc.ne.0)segact ipiloc
  284. DO 100 I=1,NIND1
  285. TYPIND = MTABTI(I)
  286. IF (TYPIND.EQ.'MOT ') THEN
  287. IP = MTABII(I)
  288. ID = IPCHAR(IP)
  289. IFI = IPCHAR(IP+1)
  290. IL1 = IFI - ID
  291. IF (IL1.EQ.ILC1) THEN
  292. IF (CHAI1.EQ.ICHARA(ID:IFI-1)) THEN
  293. IPOLY = 1
  294. ENDIF
  295. ENDIF
  296. ENDIF
  297. 100 CONTINUE
  298. if(nbesc.ne.0)segdes ipiloc
  299. SEGDES,MTABLE
  300. IF (IIMPI.EQ.333) THEN
  301. WRITE(IOIMP,*) 'DEVINI : IPOLY = ',IPOLY
  302. ENDIF
  303. IF (IPOLY.NE.0) THEN
  304. CALL ACCTAB(ITABL,'MOT',I0,X0,'VARIABLES_LIAISON_A',
  305. & L0,IP0,'TABLE',I1,X1,' ',L1,ITREFR)
  306. CALL DYNA14(ITREFR,KTLIAA)
  307. IF (IERR.NE.0) RETURN
  308. ENDIF
  309. ENDIF
  310. IF (NLIAB.NE.0) THEN
  311. CALL DEVRCO(Q1,NA1,XPTB,NPLB,XPHILB,NSB,NPLSB,NA2,IDIMB,
  312. & IBASB,IPLSB,INMSB,IORSB,3,IAROTA)
  313. CALL ACCTAB(ITABL,'MOT',I0,X0,'VARIABLES_LIAISON_B',
  314. & L0,IP0,'TABLE',I1,X1,' ',L1,ITREFR)
  315. IF (IERR.NE.0) RETURN
  316. CALL DYNE14(ITREFR,KTLIAB)
  317. IF (IERR.NE.0) RETURN
  318. ENDIF
  319. *
  320. ELSE
  321. *
  322. * Chpoints initiaux de deplacement et de vitesse:
  323. *
  324. IF (ITINIT.NE.0) THEN
  325. if (lmodyn) then
  326. mwinit = itinit
  327. segact mwinit
  328. idepl = jpdep
  329. ivite = jpvit
  330. else
  331. TYPRET=' '
  332. CALL ACCTAB(ITINIT,'MOT',I0,X0,'DEPLACEMENT',L0,IP0,
  333. & TYPRET,I1,X1,CHARRE,L1,IDEPL)
  334. *
  335. TYPRET=' '
  336. CALL ACCTAB(ITINIT,'MOT',I0,X0,'VITESSE',L0,IP0,
  337. & TYPRET,I1,X1,CHARRE,L1,IVITE)
  338. endif
  339. ENDIF
  340. *
  341. * Mise des valeurs initiales au pas m (indice 3)
  342. *
  343. IF (IDEPL.NE.0) THEN
  344. CALL DYNE18(IDEPL,KTQ,1,3,KCPR)
  345. MTQ = KTQ
  346. ENDIF
  347. IVINIT = 0
  348. IF (IVITE.NE.0) THEN
  349. CALL DYNE18(IVITE,KTQ,2,3,KCPR)
  350. IVINIT = 1
  351. ENDIF
  352. *
  353. * Ajout des forces de liaison
  354. *
  355. IF (NLIAA.NE.0) THEN
  356. CALL DEVLFA(Q1,Q2,FTOTA,NA1,IPALA,IPLIA,XPALA,XVALA,
  357. & NLIAA,XDT,1,3,FINERT,IVINIT,FTEST,FTOTA0)
  358. ENDIF
  359. IF (NLIAB.NE.0) THEN
  360. IF (RIGIDE) THEN
  361. DO 13 IP=1,NPLB
  362. DO 15 ID=1,IDIM
  363. FEXB(IP,1,ID) = FTEXB(IP,1,1,ID)
  364. 15 CONTINUE
  365. 13 CONTINUE
  366. ENDIF
  367. c SEGMENT,MTPHI
  368. c INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  369. c INTEGER IAROTA(NSB)
  370. c REAL*8 XPHILB(NSB,NPLSB,NA2,IDIMB)
  371. c ENDSEGMENT
  372. CALL DEVLF2(Q1,Q2,FTOTA,NA1,IPALB,IPLIB,XPALB,XVALB,
  373. & NLIAB,XPHILB,JPLIB,NPLB,IDIMB,FTOTB,FTOTBA,XPTB,
  374. & XDT,1,IBASB,IPLSB,INMSB,IORSB,NSB,NPLSB,NA2,3,
  375. & FEXPSM,NPC1,IERRD,FTEST2,FTOTB0,XABSCI,XORDON,
  376. & NIP,IAROTA,RIGIDE,FEXB,XCHPFB)
  377. IF (IERRD.NE.0) THEN
  378. IF (IERRD.EQ.1) CALL ERREUR(528)
  379. RETURN
  380. ENDIF
  381. ENDIF
  382. IF (NLIAA.NE.0 .OR. NLIAB.NE.0) THEN
  383. DO 5 I=1,NA1
  384. FINERT(I,4) = FINERT(I,3)
  385. 5 CONTINUE
  386. * end do
  387. ENDIF
  388. *
  389. * Calcul des deplacements et des vitesses au pas -1/2 (indice 4)
  390. *
  391. PDT = XDT(1)
  392. PDTS2 = PDT / 2.D0
  393. AUX2 = PDT * PDT / 8.D0
  394. IF (IIMPI.EQ.333) THEN
  395. DO 91 I=1,NA1
  396. WRITE(IOIMP,*)'DEVINI :'
  397. WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',3) =',FTOTA(I,3)
  398. WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',4) =',FTOTA(I,4)
  399. WRITE(IOIMP,*)'DEVINI : Fexa (',I,',1) =',Fexa(I,1)
  400. WRITE(IOIMP,*)'DEVINI : Fexa (',I,',2) =',Fexa(I,2)
  401. 91 CONTINUE
  402. * end do
  403. ENDIF
  404.  
  405. c Cas selon nature (pleine ou diag) des matrices K,C et M :
  406.  
  407. c------- Cas K,C ou M pleine -----------------------------------
  408. IF(NB1K.GT.1.OR.NB1C.GT.1.OR.NB1M.GT.1) THEN
  409. IF (IIMPI.EQ.333) WRITE(IOIMP,*)'pleine K,C,M',NB1K,NB1C,NB1M
  410.  
  411. c F_0 = F^liaison_0 + F^ext_0 - K * Q_0
  412. DO 61 I = 1,NA1
  413. FTOTA(I,3) = FTOTA(I,3) + FEXA(I,1,1)
  414. 61 CONTINUE
  415. CALL DEVLK0(Q1,XK,FTOTA,NA1,NB1K,3)
  416.  
  417. c F^AMOR_0 = C * \dot Q_0
  418. SEGINI,MAMOR
  419. CALL DEVLC0(Q2,XASM,FAMOR,NA1,NB1C,3)
  420.  
  421. c \ddot Q_0 = M-1 * (F_0 - F^AMOR_0)
  422. CALL DEVLM0(Q3,XM,XOPER,FTOTA,FAMOR,NA1,NB1,3,NB1M,3)
  423.  
  424. c Q_-1/2 = Q_0 - dt/2 * \dot Q_0 + dt^2/8 * \ddot Q_0
  425. DO 62 I = 1,NA1
  426. Q1(I,4) = Q1(I,3) - PDTS2*Q2(I,3) + AUX2*Q3(I,3)
  427. 62 CONTINUE
  428.  
  429. c F_-1/2 = F^liaison_-1/2 + F^ext_-1/2 - K * Q_-1/2
  430. DO 63 I = 1,NA1
  431. FTOTA(I,4) = FTOTA(I,4) + FEXA(I,1,2)
  432. 63 CONTINUE
  433. CALL DEVLK0(Q1,XK,FTOTA,NA1,NB1K,4)
  434.  
  435. c \dot Q_-1/2 = [4I+dt*M-1*C]-1
  436. c * ( 4 \dot Q_0 - dt*M-1*(F_-1/2+F_0-F^AMOR_0) )
  437. DO 64 I = 1,NA1
  438. FTEMP(I) = FTOTA(I,4) + FTOTA(I,3) - FAMOR(I,3)
  439. 64 CONTINUE
  440. CALL DEVLM1(QTEMP,XM,XOPER,FTEMP,NA1,NB1,3,NB1M)
  441. DO 65 I = 1,NA1
  442. FTEMP(I) = 4.D0*Q2(I,3) - PDT*QTEMP(I)
  443. 65 CONTINUE
  444. IF(NB1C.GT.1.OR.NB1M.GT.1) THEN
  445. * initialisation de l'operateur (local) et de son inverse
  446. * XOP = 4I - dt*M^-1*C
  447. SEGINI,MOP
  448. c ...C pleine, M pleine
  449. IF(NB1C.GT.1.AND.NB1M.GT.1) THEN
  450. DO 662 IA=1,NA1
  451. XOP(IA,IA) = 4.D0
  452. DO 662 IB=1,NA1
  453. AUX = 0.D0
  454. DO 663 IC=1,NA1
  455. AUX = AUX + XOPER(IA,IC,3)*XASM(IC,IB)
  456. 663 CONTINUE
  457. XOP(IA,IB) = XOP(IA,IB) - PDT*AUX
  458. 662 CONTINUE
  459. c ...C pleine, M diago
  460. ELSEIF(NB1C.GT.1) THEN
  461. DO 664 IA=1,NA1
  462. XOP(IA,IA) = 4.D0
  463. DO 664 IB=1,NA1
  464. XOP(IA,IB) = XOP(IA,IB) - PDT*XASM(IA,IB)/XM(IA,1)
  465. 664 CONTINUE
  466. c ...C diago, M pleine
  467. ELSEIF(NB1M.GT.1) THEN
  468. DO 665 IA=1,NA1
  469. XOP(IA,IA) = 4.D0
  470. DO 665 IB=1,NA1
  471. XOP(IA,IB) = XOP(IA,IB) - PDT*XOPM1(IA,IB)*XASM(IB,1)
  472. 665 CONTINUE
  473. ENDIF
  474. CALL IVMAT(NA1,XOP,INVOP,XOPM1,DETOP,0,IRET)
  475. IF (IIMPI.EQ.333) THEN
  476. WRITE(IOIMP,*) 'FTEMP(:)=',(FTEMP(jou),jou=1,NA1)
  477. do iou=1,NB1
  478. WRITE(IOIMP,*) 'XOPM1(',iou,',:)=',
  479. & (XOPM1(iou,jou),jou=1,NB1)
  480. enddo
  481. ENDIF
  482. c \dot Q_-1/2 = [4I-dt*M-1*C]-1 * FTEMP
  483. DO 67 IA = 1,NA1
  484. Q2(IA,4) = 0.D0
  485. DO 67 IB = 1,NA1
  486. Q2(IA,4) = Q2(IA,4) + XOPM1(IA,IB)*FTEMP(IB)
  487. 67 CONTINUE
  488. SEGSUP,MOP
  489. ELSE
  490. DO 68 IA = 1,NA1
  491. Q2(IA,4) = FTEMP(IA) / (4.D0 - (PDT*XASM(IA,1)/XM(IA,1)))
  492. 68 CONTINUE
  493. ENDIF
  494.  
  495. c \ddot Q_-1/2 = M-1 * ( F_-1/2 - F^AMOR_-1/2)
  496. CALL DEVLC0(Q2,XASM,FAMOR,NA1,NB1C,4)
  497. CALL DEVLM0(Q3,XM,XOPER,FTOTA,FAMOR,NA1,NB1,3,NB1M,4)
  498.  
  499. SEGSUP,MAMOR
  500.  
  501. c------- Cas K,C et M diagonal -----------------------------------
  502. ELSE
  503.  
  504. DO 10 I = 1,NA1
  505. UNSM3 = 1.D0 / ( XM(I,1) - FINERT(I,3) )
  506. UNSM4 = 1.D0 / ( XM(I,1) - FINERT(I,4) )
  507. *
  508. * G = F - K Q
  509. * i,0 i,0 i i,0
  510. *
  511. FTOTA(I,3) = FTOTA(I,3) + FEXA(I,1,1) - (XK(I,1)*Q1(I,3))
  512. AUX1 = (FTOTA(I,3)*UNSM3) - (XASM(I,1)*UNSM3*Q2(I,3))
  513. AUX3 = Q1(I,3) - ( PDTS2 * Q2(I,3) )
  514. *
  515. * . 2 .
  516. * Q = Q - h/2 Q + h /8 ( G - A Q )
  517. * i,-1/2 i,0 i,0 i,0 i i,0
  518. *
  519. Q1(I,4) = AUX3 + ( AUX2 * AUX1 )
  520. *
  521. * G = F - K Q
  522. * i,-1/2 i,-1/2 i i,-1/2
  523. *
  524. FTOTA(I,4) = FTOTA(I,4) + FEXA(I,1,2) - (XK(I,1)*Q1(I,4))
  525. AUX4 = 4.D0 + ( XASM(I,1) * PDT * UNSM3 )
  526. AUX5 = ( FTOTA(I,4) * UNSM4 ) + ( FTOTA(I,3) * UNSM3 )
  527. AUX6 = 1.D0 / ( 4.D0 - ( XASM(I,1) * PDT * UNSM3 ) )
  528. *
  529. * . 1 .
  530. * Q = ________ ( ( 4 + A h ) Q - h ( G + G ) )
  531. * i,-1/2 4 - A h i i,0 i,-1/2 i,0
  532. * i
  533. *
  534. Q2(I,4) = AUX6 * ( (AUX4*Q2(I,3)) - (PDT*AUX5) )
  535. *
  536. * " .
  537. * Q = ( G - A Q ) / M
  538. * i i i i i
  539. *
  540. Q3(I,3) = (FTOTA(I,3)*UNSM3) - (XASM(I,1)*UNSM3*Q2(I,3))
  541. Q3(I,4) = (FTOTA(I,4)*UNSM4) - (XASM(I,1)*UNSM4*Q2(I,4))
  542. 10 CONTINUE
  543.  
  544. ENDIF
  545. *
  546. ENDIF
  547. IF (IIMPI.EQ.333) THEN
  548. DO 30 I=1,NA1
  549. WRITE(IOIMP,*)'DEVINI :'
  550. WRITE(IOIMP,*)'DEVINI : Q1(',I,',3) =',Q1(I,3)
  551. WRITE(IOIMP,*)'DEVINI : Q2(',I,',3) =',Q2(I,3)
  552. WRITE(IOIMP,*)'DEVINI : Q3(',I,',3) =',Q3(I,3)
  553. WRITE(IOIMP,*)'DEVINI :'
  554. WRITE(IOIMP,*)'DEVINI : Q1(',I,',4) =',Q1(I,4)
  555. WRITE(IOIMP,*)'DEVINI : Q2(',I,',4) =',Q2(I,4)
  556. WRITE(IOIMP,*)'DEVINI : Q3(',I,',4) =',Q3(I,4)
  557. 30 CONTINUE
  558. * end do
  559. ENDIF
  560. IF (IIMPI.EQ.333) THEN
  561. DO 31 I=1,NA1
  562. WRITE(IOIMP,*)'DEVINI :'
  563. WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',3) =',FTOTA(I,3)
  564. WRITE(IOIMP,*)'DEVINI : FTOTA(',I,',4) =',FTOTA(I,4)
  565. 31 CONTINUE
  566. * end do
  567. ENDIF
  568.  
  569. *
  570. ICPR = KCPR
  571. SEGSUP,ICPR
  572. *
  573. END
  574.  
  575.  
  576.  
  577.  
  578.  
  579.  
  580.  
  581.  
  582.  

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