Télécharger mondes.eso

Retour à la liste

Numérotation des lignes :

mondes
  1. C MONDES SOURCE PV090527 24/01/19 21:15:05 11827
  2. SUBROUTINE MONDES(MMATRX,MVECTX,NOEN,ISOUCI)
  3. C
  4. C **** EXECUTE LA SOLUTION X DE (Lt D L) X=F
  5. C
  6. CMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMB
  7. CMB
  8. CMB Plutot la solution de L.D.Lt ou L.D.Mt (cas non symétrique)
  9. CMB Elle devrait dons s'appeller DESMON et non MONDES.
  10. CMB
  11. CMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMBCMB
  12. C
  13. IMPLICIT INTEGER(I-N)
  14. IMPLICIT REAL*8(A-H,O-Z)
  15. PARAMETER (LPCRAY=10000)
  16. INTEGER OOOVAL
  17.  
  18. C ... Variable décrivant l'EXIStence du Triangle Supérieur différent
  19. C de l'inférieur (cas non symétrique) ...
  20. LOGICAL EXISTS
  21.  
  22. -INC SMMATRI
  23. -INC SMELEME
  24. -INC SMVECTD
  25.  
  26. -INC PPARAM
  27. -INC CCOPTIO
  28. -INC CCREEL
  29. -INC SILICRE
  30. -INC CCASSIS
  31.  
  32. SEGMENT ITTRV(MULRE)
  33. SEGMENT,ITRAA(NENSL)
  34. SEGMENT/BID/(IBIDON(IIMAX))
  35. logical toutmem
  36. toutmem=.true.
  37. DNORMA=0.D0
  38. DNORMB=0.D0
  39. crit=xpetit
  40. crit2 = 0.D0
  41. call oooprl(1)
  42. nbthr=nbthrs
  43. if (noen.eq.0) then
  44. nbelem=0
  45. nbnn =1
  46. nbsous=0
  47. nbref =0
  48. segini,ipt8
  49. ipt8.itypel=1
  50. endif
  51. CALL OOOMRU(1)
  52. IVALMA=0
  53. IF(IIMPI.EQ.3) WRITE(IOIMP,1000)MMATRX,MVECTX
  54. 1000 FORMAT(' SUBROUTINE MONDES : POINTEUR DE LA MATRICE=',I5,
  55. 1 ' POINTEUR DU VECTEUR=',I5)
  56. C
  57. C **** ACTIVATION DES SEGMENTS
  58. C
  59. MMATRI=MMATRX
  60. SEGACT,MMATRI*MOD
  61. * BEC00SE OPTIMISEUR
  62. ITRAA=MMATRI
  63. NENSL=0
  64. IF(NENS.NE.0) THEN
  65. NENSL=NENS
  66. SEGINI,ITRAA
  67. ENDIF
  68. NE1=NENSL
  69.  
  70. EXISTS=.FALSE.
  71. IF(IILIGS.NE.0) THEN
  72. MILIGN=IILIGS
  73. SEGACT MILIGN
  74. IF(ILIGN(/1).GT.0) EXISTS=.TRUE.
  75. ENDIF
  76.  
  77. MILIGN=IILIGN
  78. SEGACT,MILIGN
  79. INC=IPNO(/1)
  80. nbthr=min(nbthrs,inc/1000+1)
  81.  
  82. MVECTD=MVECTX
  83. SEGACT MVECTD*MOD
  84.  
  85. MDNOR=IDNORM
  86. SEGACT MDNOR
  87. IF(IDNORD.GT.0) THEN
  88. MDNO1=IDNORD
  89. SEGACT MDNO1
  90. ELSE
  91. MDNO1=MDNOR
  92. ENDIF
  93.  
  94. MDIAG=IDIAG
  95.  
  96. C ... MULRE = nombre de seconds membres ...
  97. MULRE = VECTBB(/1)/INC
  98. C ... MUINC ne servira que comme une borne sur les indices de boucles ...
  99. MUINC = ( MULRE-1)*INC
  100. SEGINI ITTRV
  101. inc=vectbb(/1)
  102. * mvect1 sauve les forces
  103. * mvect2 sauve le deplacement
  104. * mvect3 vecteur de travail
  105. * mvectd evolue: force a resoudre (residu) puis increment de deplacement
  106. if (mulre.eq.1.and..not.exists) segini mvect1,mvect2
  107. INC=IPNO(/1)
  108.  
  109. C ... Multiplication du second membre par les DNOR(.) (à cause du
  110. C conditionnement des matrices de blocages) ...
  111. DNORMA=0.D0
  112. DO 450 K=0,MUINC,INC
  113. DO 45 I=1,INC
  114. VECTBB(I+K)=VECTBB(I+K)*MDNO1.DNOR(I)
  115. ** DNORMA=MAX(DNORMA,ABS(VECTBB(I+K)))
  116. 45 CONTINUE
  117. 450 CONTINUE
  118. ** DNORMA=DNORMA*xzprec*xzprec
  119. DNORMA=xpetit/xzprec
  120.  
  121. C ... Détection du premier terme important du second membre ...
  122. II=0
  123. DO 451 K=0,MUINC,INC
  124. II=II+1
  125. DO 452 I=1,INC
  126. IF(ABS(VECTBB(I+K)).LT.DNORMA) GO TO 452
  127. C ... Le numéro du noeud concerné par le premier terme important va à ITTRV ...
  128. ITTRV(II)=IPNO(I)
  129. GO TO 451
  130. 452 CONTINUE
  131. 451 CONTINUE
  132.  
  133. C ... NNOE = nombre de noeuds concernés ...
  134. NNOE=ILIGN(/1)
  135. IDEP=NNOE
  136. C ... On cherche le minimum (non nul) des ITTRV que l'on met dans IDEP ...
  137. DO 453 I=1,MULRE
  138. IF(ITTRV(I).LT.IDEP) THEN
  139. IF(ITTRV(I).NE.0) THEN
  140. IDEP=ITTRV(I)
  141. ELSE
  142. IDEP=1
  143. GO TO 4530
  144. ENDIF
  145. ENDIF
  146. 453 CONTINUE
  147. 4530 CONTINUE
  148. SEGDES MDNOR
  149. *
  150. * debut de la boucle d'amelioration de la precision si necessaire
  151. *
  152. * sauver le champs de forces initial
  153. if(mulre.eq.1.and..not.exists) then
  154. do i=1,inc
  155. mvect1.vectbb(i)=vectbb(i)
  156. enddo
  157. icorrec=0
  158. endif
  159. 10000 continue
  160. IAA=0
  161.  
  162.  
  163.  
  164.  
  165.  
  166. SEGACT,MDIAG
  167. C
  168. C **** DESCENTE: ON RESOU L*C=B. EN FAIT ON STOCKE C DANS B
  169. C
  170. LTRK=MAX(1+LPCRAY,OOOVAL(1,4))
  171. IIMAX=(((IJMAX +LTRK)/LTRK)+1)*LTRK
  172. iimax=oooval(1,1)/10
  173. * test si la place disponible permet de tout mettre en memoire
  174. xplds=oooval(1,1)-oooval(2,3)
  175. if (real(nnoe)*ijmax.lt.xplds) iimax=0
  176. SEGINI /err=3500/ BID
  177. goto 3501
  178. 3500 continue
  179. * place insuffisante en memoire centrale. On en fait
  180. * desactivation de la matrice en partant du debut
  181. do ivs=1,nnoe
  182. lign=ilign(ivs)
  183. segdes lign
  184. enddo
  185. segini bid
  186. 3501 continue
  187. C ... IDEP : on commence par ce noeud car le second membre est nul pour
  188. C tous les precedents ...
  189. * activation de tous les derniers noeuds car ils sont deja en memoire normalement
  190. do 2001 ivs=nnoe,1,-1
  191. lign=ilign(ivs)
  192. segact/err=2010/ lign
  193. 2001 continue
  194. 2010 continue
  195. ilmin=ivs+1
  196. DO 610 IVS=IDEP,NNOE
  197. LIGN=ILIGN(IVS)
  198. SEGACT /ERR=50/ LIGN
  199. IVALMA=IVALMA+VAL(/1)
  200. *pvpv IF (IVALMA.GT.NGMAXY) GOTO 50
  201. NA=IMMM(/1)
  202. IPRELL=IPREL
  203. DO 611 IA=1,NA
  204. C ... Si l'inconnue présente un mode d'ensemble ...
  205. IF(IMMM(IA).NE.0) THEN
  206. C ... On incrémente le compteur et
  207. IAA=IAA+1
  208. C ... on met le N° de la ligne dans ITRAA ...
  209. ITRAA(IAA)=IPRELL+IA-1
  210. ENDIF
  211. 611 CONTINUE
  212. segact lign
  213. CALL MONDE2(ITTRV(1),IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  214. > NA,IPREL,MULRE,INC,IVS,LLOM,dnorma)
  215. 610 CONTINUE
  216. C ... Si on n'a pas eu de pb de mémoire on passe par là ...
  217. SEGSUP BID
  218. C ... ILMAX vaut le dernier noeud qui a pu être traité ...
  219. ILMAX=NNOE
  220. C ... On va à la division par la diagonale ...
  221. GOTO 55
  222.  
  223. 50 CONTINUE
  224. C ... Si on est là, c'est à cause des pb de mémoire ...
  225. SEGSUP BID
  226. toutmem=.false.
  227. C ... IVS est le N° du LIGN qui n'a pas pu être traité ...
  228. C ... ILMAX = N° du dernier traité ...
  229. ILMAX=IVS-1
  230. C ... IILMAX = N° du premier à traiter ...
  231. IILMAX=IVS
  232. DO 54 IVS=IILMAX,NNOE
  233. 58 CONTINUE
  234. LIGN=ILIGN(IVS)
  235. C ... Même si tout à l'heure SEGACT n'a pas marché, maintenant on a supprimé BID ...
  236. SEGACT /ERR=56/ LIGN
  237. GOTO 57
  238. 56 CONTINUE
  239. * EN CAS DE PROBLEME ON FAIT UN PEU DE PLACE
  240. toutmem=.false.
  241. C ... Si on a toujours des pb de mémoire, on désactive le LIGN N° ILMAX,
  242. LIGN=ILIGN(ILMAX)
  243. if(ilmax.lt.ilmin) SEGDES LIGN*(NOMOD,MRU)
  244. C ... puis ce dernier est décrémenté ...
  245. ILMAX=ILMAX-1
  246. IF (ILMAX.EQ.1) CALL ERREUR(5)
  247. C ... et on recommence ...
  248. GOTO 58
  249. C ... On est là si tout s'est bien passé avec SEGACT ...
  250. 57 CONTINUE
  251. NA=IMMM(/1)
  252. IPRELL=IPREL
  253. DO 612 IA=1,NA
  254. IF(IMMM(IA).EQ.0) GOTO 612
  255. IAA=IAA+1
  256. ITRAA(IAA)=IPRELL+IA-1
  257. 612 CONTINUE
  258. CALL MONDE2(ITTRV(1),IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  259. > NA,IPREL,MULRE,INC,IVS,LLOM,dnorma)
  260. if (ivs.lt.ilmin) SEGDES,LIGN*(NOMOD,MRU)
  261. 54 CONTINUE
  262. 55 CONTINUE
  263. C
  264. C ... À cet endroit la descente est finie. L'état des LIGN est le suivant :
  265. C ... Ceux de IDEP à ILMAX sont actifs, les autres (si ILMAX < NNOE) sont désactivés ...
  266. C
  267. C **** DIVISION PAR LE TERME DIAGONAL ****
  268. C
  269. dnorma=0.d0
  270. idenorm=0
  271. DO 120 K=0,MUINC,INC
  272. DO 12 I=1,INC
  273. J=I+K
  274. if (2*abs(diag(i)).gt.abs(diag(i))) goto 122
  275. idenorm=1
  276. 122 continue
  277. VECTBB(J)=VECTBB(J)*DIAG(I)
  278. ** dnorma=max(dnorma,abs(vectbb(j)))
  279. 12 CONTINUE
  280.  
  281. 120 CONTINUE
  282. ** dnorma=dnorma*xzprec*xzprec
  283. dnorma=xpetit/xzprec
  284. SEGDES,MDIAG
  285.  
  286. C
  287. C **** MONTEE ****
  288. C
  289. C ... Si la matrice est asymétrique ...
  290. IF(EXISTS) THEN
  291.  
  292. C ... On commence par désactiver les LIGN qui sont restés actifs ...
  293. DO 2000 I=IDEP,ILMAX
  294. LIGN=ILIGN(I)
  295. if (i.lt.ilmin) seGDES,LIGN*(NOMOD,MRU)
  296. 2000 CONTINUE
  297.  
  298. C ... Puis on désactive MILIGN ...
  299. SEGDES,MILIGN
  300.  
  301. C ... On change de MILIGN ...
  302. MILIGN=IILIGS
  303. SEGACT,MILIGN
  304.  
  305. C ... Et on active des LIGN pointés par ce dernier ...
  306. ILMAX=IDEP-1
  307. DO 2100 I=IDEP,NNOE
  308. LIGN=ILIGN(I)
  309. SEGACT /ERR=2110/ LIGN
  310. ILMAX = I
  311. 2100 CONTINUE
  312. C ... Et on passe à la montée ...
  313. GOTO 3000
  314.  
  315. C ... Si on n'a même pas activé le premier, alors ADIOS ...
  316. 2110 IF(ILMAX.EQ.IDEP-1) CALL ERREUR(5)
  317. toutmem=.false.
  318. 3000 CONTINUE
  319. ENDIF
  320.  
  321.  
  322. J=NNOE
  323. C ... Début de la pseudo boucle sur les LIGN qui sont restés désactivés ...
  324. 70 CONTINUE
  325. IF (J.EQ.ILMAX) GOTO 72
  326. LIGN=ILIGN(J)
  327. SEGACT /ERR=68/ LIGN
  328. GO TO 67
  329. C ... Si on a des pb avec activation, on désactive et on diminue ILMAX ...
  330. 68 CONTINUE
  331. toutmem=.false.
  332. LIGN = ILIGN( ILMAX)
  333. if(ilmax.lt.ilmin) SEGDES LIGN*(NOMOD,MRU)
  334. ILMAX=ILMAX-1
  335. IF(ILMAX.EQ.1) CALL ERREUR (5)
  336. C ... puis on recommence la tentative ...
  337. GO TO 70
  338. C ... On a réussi à activer ...
  339. 67 CONTINUE
  340. NA=IMMM(/1)
  341. IPRELL=IPREL
  342. CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  343. > NA,IPREL,MULRE,INC,dnorma)
  344. if (j.lt.ilmin) SEGDES,LIGN*(NOMOD,MRU)
  345. J = J-1
  346. GO TO 70
  347. 72 CONTINUE
  348. CC FIN DE PSEUDO BOUCLE J = INC ,ILMAX+1,-1 - Vieux commentaire et Faux !!!
  349. CC FIN DE PSEUDO BOUCLE J = NNOE ,ILMAX+1,-1
  350.  
  351. C ... Dans cette boucle on parcourt les LIGN qui sont restés actifs ...
  352. DO 13 J=ILMAX,IDEP,-1
  353. LIGN=ILIGN(J)
  354. NA=IMMM(/1)
  355. IPRELL=IPREL
  356. CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  357. > NA,IPREL,MULRE,INC,dnorma)
  358. if (j.lt.ilmin) SEGDES,LIGN*(NOMOD,MRU)
  359. 13 CONTINUE
  360.  
  361. C ... On n'avait même pas essayé d'activer les IDEP-1 premiers LIGN ...
  362. DO 113 J=IDEP-1,1,-1
  363. LIGN=ILIGN(J)
  364. SEGACT LIGN
  365. NA=IMMM(/1)
  366. IPRELL=IPREL
  367. DO 1131 ILM=1,NA
  368. IF(IMMM(ILM).EQ.0) GO TO 1131
  369. IAA=IAA+1
  370. ITRAA(IAA)=IPRELL+ILM-1
  371. 1131 CONTINUE
  372. CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  373. > NA,IPREL,MULRE,INC,dnorma)
  374. if (j.lt.ilmin) SEGDES,LIGN*(NOMOD,MRU)
  375. 113 CONTINUE
  376.  
  377. if (idenorm.eq.1) then
  378. call erreur(1049)
  379. do k=0,muinc,inc
  380. do i=1,inc
  381. vectbb(i+k)=0.d0
  382. enddo
  383. enddo
  384. endif
  385. * desactivation generale
  386. do ivs=ilmin,nnoe
  387. lign=ilign(ivs)
  388. segdes lign
  389. enddo
  390.  
  391. *
  392. * on verifie la precision du resultat si pas de noen ou pas de modes d'ensemble
  393. ** write(6,*) ' noen nens ',noen,nens
  394. if (noen.ne.0.or.nens.eq.0) then
  395. if(mulre.eq.1.and..not.exists) then
  396. * U+dU
  397. if (icorrec.eq.0) then
  398. do i=1,inc
  399. mvect2.vectbb(i)=vectbb(i)
  400. enddo
  401. else
  402. do i=1,inc
  403. mvect2.vectbb(i)=mvect2.vectbb(i)+vectbb(i)
  404. enddo
  405. endif
  406. ilicre=jlicre
  407. segact ilicre
  408. ligcre=ligcrp
  409. segact ligcre
  410. * K*U
  411. segini mvect3
  412. call graco7(ilicre,mvect2,mvect3,inc,nbthr,0)
  413. * F-K*U
  414. crit=xpetit/xzprec
  415. crit2 = 0.D0
  416. do i=1,inc
  417. crit=max(crit,abs(mvect1.vectbb(i)),abs(mvect3.vectbb(i)))
  418. mvect3.vectbb(i)=mvect1.vectbb(i)-mvect3.vectbb(i)
  419. ** write(6,*) ' mvect1 mvect3 ',mvect1.vectbb(i),mvect3.vectbb(i)
  420. crit2=max(crit2,abs(mvect3.vectbb(i)))
  421. enddo
  422. * test
  423. if (icorrec.le.3) crite = crit * sqrt(xzprec*xszpre)
  424.  
  425. * critere lache apres les premieres ite©rations. demande une reflexion approfondie.
  426. if (icorrec.ge.5) crite = crit * 1d-4
  427. if(icorrec.gt.7.and.nensl.ne.0) crite=xgrand*xzprec
  428. ** write(6,*) ' crite crit2',crite,crit2,inc
  429. if (crit2.gt.crite) then
  430. * il faut iterer
  431. do i2=1,inc
  432. vectbb(i2)=mvect3.vectbb(i2)/(icorrec+1)
  433. enddo
  434. segsup mvect3
  435. icorrec=icorrec+1
  436. if (icorrec.lt.10.and.toutmem) goto 10000
  437. endif
  438. * verif correcte
  439. segsup mvect1,mvect3
  440. do i=1,inc
  441. vectbb(i)=mvect2.vectbb(i)
  442. enddo
  443. segsup mvect3
  444. reaerr(1)=crit2/crit
  445. interr(1)=icorrec
  446. interr(2)=nensl
  447. if (icorrec.ge.10) then
  448. if (isouci.eq.0) then
  449. call erreur(1128)
  450. else
  451. call soucis(1128)
  452. endif
  453. endif
  454. if (icorrec.ge.1) then
  455. * l'erreur 1129 est informative donc pas de souci
  456. call erreur(1129)
  457. endif
  458. endif
  459. endif
  460. ** write(6,*) 'icorrec',icorrec
  461. C
  462. C **** ON VERIFIE QUE PAS DE MODE RIGIDE ACTIVE
  463. C
  464. C ... IAA = nombre de modes d'ensemble (selon les indications dans IMMM) ...
  465. NENS=IAA
  466. ** write(6,*) ' mondes nens nensl ',nens,nensl
  467. NBEN=0
  468. IF(NENS.EQ.0) GO TO 26
  469.  
  470. MINCPO=IINCPO
  471. MIMIK=IIMIK
  472. SEGACT,MINCPO,MIMIK
  473. MELEME=IGEOMA
  474. SEGACT,MELEME
  475. NNOE=INCPO(/2)
  476. IINC1=INCPO(/1)
  477.  
  478. C ... Boucle sur les seconds membres ...
  479. DO 300 KI=0,MUINC,INC
  480. C ... XMA et XMAL seront le maximum des val. abs. des termes
  481. C du résultat partiel (avant la multiplication par MDNOR)
  482. C N° KI+1 multiplié par 1.e-10 ...
  483. XMA=xpetit/xzprec
  484. XMAL=xpetit/xzprec
  485. DO 30 KK=1,INC
  486. I=KK+KI
  487. if (ittr(kk).eq.0) then
  488. XMA=MAX(XMA,ABS(VECTBB(I)))
  489. else
  490. XMAL=MAX(XMAL,ABS(VECTBB(I)))
  491. endif
  492. 30 CONTINUE
  493. XMA=XMA*1.d-10
  494. XMAL=XMAL*max(1d-10,crit2/crit*1d5)
  495. xmal=max(xma*1d-2,xmal)
  496. * write (6,*) ' mondes xma xmal ',xma,xmal
  497.  
  498. C ... Boucle sur les modes d'ensemble ...
  499. ** write (6,*) ' nombre de modes d ensemble',nens
  500. iwrite = 0
  501. DO 21 IA=1,NENS
  502. I1=ITRAA(IA)
  503. J=IPNO(I1)
  504. DO 22 K=1,IINC1
  505. IF(INCPO(K,J).EQ.I1) GO TO 23
  506. 22 CONTINUE
  507. CALL ERREUR(5)
  508. call oooprl(0)
  509. RETURN
  510.  
  511. 23 CONTINUE
  512. C ... Si ce n'est pas un multiplicateur, le déplacement doit être
  513. C < à XMA, sinon ERREUR ...
  514. * write (6,*) ' mondes i1 ittr vect xma ',
  515. * > i1,ittr(i1),vectbb(i1+ki),xma
  516. IF(ITTR(I1).EQ.0.AND.ABS(VECTBB(I1+KI)).LE.XMA ) then
  517. vectbb(i1+ki)=0.d0
  518. GO TO 20
  519. endif
  520. C ... Si c'est un multiplicateur, le multiplicateur doit être
  521. C < à XMAL sinon ERREUR ...
  522. if (ittr(i1).ne.0) then
  523. i2=ittr(i1)
  524. vecs=vectbb(i1+ki)+vectbb(i2+ki)
  525. if (abs(vecs).le.xmal) then
  526. vectbb(i1+ki)=0.d0
  527. vectbb(i2+ki)=0.d0
  528. goto 20
  529. endif
  530. endif
  531. ** write (6,*) ' ittr vect ',i1,i2,vectbb(i1+ki),
  532. ** > vectbb(i2+ki),xmal
  533.  
  534.  
  535. C ... Si option NOEN on ne fait pas d'erreur, on accumule les pts dans un meleme
  536. C ... Si option NOEN on accumule les pts dans un meleme
  537. IF(NOEN.EQ.0) THEN
  538. nbelem=ipt8.num(/2)+1
  539. segadj ipt8
  540. ipt8.num(1,nbelem)=num(1,j)
  541. ENDIF
  542. NBEN=NBEN+1
  543. MOTERR(1:4)=IMIK(K)
  544. INTERR(1)=NUM(1,J)
  545. IF(NOEN.EQ.0.OR.ISOUCI.EQ.1) THEN
  546. call soucis(149)
  547. GO TO 21
  548. ENDIF
  549. CALL ERREUR(149)
  550. call oooprl(0)
  551. RETURN
  552.  
  553. 20 CONTINUE
  554. C ... Messages d'information ...
  555. JJK=NUM(1,J)
  556. KIK=KI/INC +1
  557.  
  558.  
  559. * on n'ecrit qu'une seule fois le message indetermination
  560. IF(ITTR(I1).EQ.0 .AND. MULRE.EQ.1.and.iwrite.eq.0) then
  561. if (imik(k).ne.'LX ')
  562. > WRITE(IOIMP,41) JJK,IMIK(K)
  563. iwrite=iwrite+1
  564. endif
  565.  
  566. IF (IIMPI.NE.0.AND. ITTR(I1).NE.0.AND.MULRE.EQ.1) then
  567. WRITE(IOIMP,42) JJK,IMIK(K)
  568. endif
  569.  
  570. IF(ITTR(I1).eq.0 .AND. MULRE.NE.1.and.iwrite.eq.0) then
  571. WRITE(IOIMP,410)KIK,JJK,IMIK(K)
  572. iwrite=iwrite+1
  573. endif
  574. IF (IIMPI.NE.0 .AND. ITTR(I1).NE.0.AND.MULRE.NE.1)
  575. & WRITE(IOIMP,420)KIK,JJK,IMIK(K)
  576.  
  577. 21 CONTINUE
  578.  
  579. 300 CONTINUE
  580.  
  581. 41 FORMAT(' INDETERMINATION DETECTEE AU NOEUD ',I6,' INCONNUE ',
  582. * A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ',
  583. * 'LA SUSDITE dans mondes')
  584. 42 FORMAT(' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ',
  585. * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ',
  586. * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE dans mondes')
  587. 410 FORMAT(' VECTEUR NUMERO',I3,' INDETERMINATION DETECTE AU NOEUD '
  588. *,I6,' INCONNUE ',
  589. * A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ',
  590. * 'LA SUSDITE dans mondes')
  591. 420 FORMAT(' VECTEUR NUMERO ',I3,/,
  592. * ' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ',
  593. * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ',
  594. * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE dans mondes')
  595.  
  596. SEGDES,MELEME
  597. SEGDES,MINCPO
  598. SEGDES,MIMIK
  599. 26 CONTINUE
  600. if (noen.eq.0) then
  601. * test si nan ou inf dans le resultat
  602. inan = 0
  603. DO 500 KI=0,MUINC,INC
  604. DO 501 KK=1,INC
  605. i = kk + ki
  606. if (abs(vectbb(i)).lt.xgrand) goto 501
  607. inan = inan + 1
  608. vectbb(i)=xgrand
  609. 501 continue
  610. 500 continue
  611. if (inan.ne.0) then
  612. nben = -inan
  613. call soucis(nben)
  614. endif
  615. endif
  616. * indiquer si besoin le nombre de modes d'ensemble excites
  617. * et le maillage des noeuds freinés
  618. IF (NOEN.EQ.0) then
  619. * write (6,*) ' mondes nben ',nben
  620. CALL ECRENT(NBEN)
  621. call ecrobj('MAILLAGE',ipt8)
  622. segdes ipt8
  623. endif
  624. MDNOR=IDNORM
  625. SEGACT,MDNOR
  626. DO 350 K=0,MUINC,INC
  627. DO 35 I = 1, INC
  628. VECTBB(I+K)=VECTBB(I+K)*DNOR(I)
  629. 35 CONTINUE
  630. * on force l'egalite des multiplicateurs de lagrange
  631. DO 36 I = 1, INC
  632. if (ITTR(I).ne.0) then
  633. * write (6,*) ' dans mondes ',i,ittr(i)
  634. if (vectbb(i+k).eq.0.d0.or.vectbb(ittr(i)+k).eq.0.d0) then
  635. * write (6,*) ' mondes vectbbs ',vectbb(i+k),vectbb(ittr(i)+k)
  636. vectbb(i+k)=0.d0
  637. vectbb(ittr(i)+k)=0.d0
  638. goto 36
  639. endif
  640. vectbb(i+k)=(vectbb(i+k)+vectbb(ittr(i)+k))/2
  641. vectbb(ittr(i)+k)=vectbb(i+k)
  642. endif
  643. 36 CONTINUE
  644. 350 CONTINUE
  645. SEGDES,MDNOR
  646. C ... On ne désactive MDNO1 que dans le cas où il est
  647. C différent de MDNOR ...
  648. IF(IDNORD.GT.0) THEN
  649. SEGDES,MDNO1
  650. ENDIF
  651. SEGDES MMATRI
  652. SEGDES,MILIGN
  653. SEGDES,MVECTD
  654. SEGSUP ITTRV
  655. IF(NENSL.NE.0) SEGSUP,ITRAA
  656.  
  657. IF(IIMPI.EQ.3) WRITE(IOIMP,1001) MVECTD
  658. 1001 FORMAT(' SUBROUTINE MONDES : POINTEUR DU VECTEUR EN SORTIE=',I5)
  659.  
  660. CALL OOOMRU(0)
  661. call oooprl(0)
  662. RETURN
  663. END
  664.  
  665.  
  666.  
  667.  
  668.  
  669.  
  670.  
  671.  
  672.  
  673.  
  674.  
  675.  
  676.  
  677.  
  678.  
  679.  
  680.  
  681.  
  682.  
  683.  
  684.  
  685.  
  686.  
  687.  
  688.  

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