Télécharger mondes.eso

Retour à la liste

Numérotation des lignes :

  1. C MONDES SOURCE CB215821 17/05/22 21:15:03 9438
  2. SUBROUTINE MONDES(MMATRX,MVECTX,NOEN)
  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. -INC CCOPTIO
  26. -INC CCREEL
  27.  
  28. SEGMENT ITTRV(MULRE)
  29. SEGMENT,ITRAA(NENS)
  30. SEGMENT/BID/(BIDON(IIMAX+10)*D)
  31.  
  32. DNORMA=0.D0
  33. DNORMB=0.D0
  34. if (noen.eq.0) then
  35. nbelem=0
  36. nbnn =1
  37. nbsous=0
  38. nbref =0
  39. segini,ipt8
  40. ipt8.itypel=1
  41. endif
  42. CALL OOOMRU(1)
  43. IVALMA=0
  44. IF(IIMPI.EQ.3) WRITE(IOIMP,1000)MMATRX,MVECTX
  45. 1000 FORMAT(' SUBROUTINE MONDES : POINTEUR DE LA MATRICE=',I5,
  46. 1 ' POINTEUR DU VECTEUR=',I5)
  47. C
  48. C **** ACTIVATION DES SEGMENTS
  49. C
  50. MMATRI=MMATRX
  51. SEGACT,MMATRI*MOD
  52. * BEC00SE OPTIMISEUR
  53. ITRAA=MMATRI
  54. IF(NENS.NE.0) SEGINI,ITRAA
  55. NE1=NENS
  56. IAA=0
  57.  
  58. EXISTS=.FALSE.
  59. IF(IILIGS.NE.0) THEN
  60. MILIGN=IILIGS
  61. SEGACT MILIGN
  62. IF(ILIGN(/1).GT.0) EXISTS=.TRUE.
  63. ENDIF
  64.  
  65. MILIGN=IILIGN
  66. SEGACT,MILIGN
  67. INC=IPNO(/1)
  68.  
  69. MVECTD=MVECTX
  70. SEGACT MVECTD*MOD
  71.  
  72. MDNOR=IDNORM
  73. SEGACT MDNOR
  74. IF(IDNORD.GT.0) THEN
  75. MDNO1=IDNORD
  76. SEGACT MDNO1
  77. ELSE
  78. MDNO1=MDNOR
  79. ENDIF
  80.  
  81. MDIAG=IDIAG
  82. SEGACT,MDIAG
  83.  
  84. C ... MULRE = nombre de seconds membres ...
  85. MULRE = VECTBB(/1)/INC
  86. C ... MUINC ne servira que comme une borne sur les indices de boucles ...
  87. MUINC = ( MULRE-1)*INC
  88. SEGINI ITTRV
  89.  
  90. C ... Multiplication du second membre par les DNOR(.) (à cause du
  91. C conditionnement des matrices de blocages) ...
  92. DNORMA=0.D0
  93. DO 450 K=0,MUINC,INC
  94. DO 45 I=1,INC
  95. VECTBB(I+K)=VECTBB(I+K)*MDNO1.DNOR(I)
  96. DNORMA=MAX(DNORMA,ABS(VECTBB(I+K)))
  97. 45 CONTINUE
  98. 450 CONTINUE
  99. DNORMA=DNORMA*xzprec*xzprec
  100.  
  101. C ... Détection du premier terme important du second membre ...
  102. II=0
  103. DO 451 K=0,MUINC,INC
  104. II=II+1
  105. DO 452 I=1,INC
  106. IF(ABS(VECTBB(I+K)).LT.DNORMA) GO TO 452
  107. C ... Le numéro du noeud concerné par le premier terme important va à ITTRV ...
  108. ITTRV(II)=IPNO(I)
  109. GO TO 451
  110. 452 CONTINUE
  111. 451 CONTINUE
  112.  
  113. C ... NNOE = nombre de noeuds concernés ...
  114. NNOE=ILIGN(/1)
  115. IDEP=NNOE
  116. C ... On cherche le minimum (non nul) des ITTRV que l'on met dans IDEP ...
  117. DO 453 I=1,MULRE
  118. IF(ITTRV(I).LT.IDEP) THEN
  119. IF(ITTRV(I).NE.0) THEN
  120. IDEP=ITTRV(I)
  121. ELSE
  122. IDEP=1
  123. GO TO 4530
  124. ENDIF
  125. ENDIF
  126. 453 CONTINUE
  127. 4530 CONTINUE
  128. SEGDES MDNOR
  129. C
  130. C **** DESCENTE: ON RESOU L*C=B. EN FAIT ON STOCKE C DANS B
  131. C
  132. LTRK=MAX(1+LPCRAY,OOOVAL(1,4))
  133. IIMAX=(((IJMAX +LTRK)/LTRK)+1)*LTRK
  134. * test si la place disponible permet de tout mettre en memoire
  135. xplds=oooval(1,1)-oooval(2,3)
  136. if (real(nnoe)*ijmax.lt.xplds) iimax=0
  137. SEGINI BID
  138. C ... IDEP : on commence par ce noeud car le second membre est nul pour
  139. C tous les précédents ...
  140. DO 610 IVS=IDEP,NNOE
  141. LIGN=ILIGN(IVS)
  142. SEGACT /ERR=50/ LIGN
  143. IVALMA=IVALMA+VAL(/1)
  144. *pvpv IF (IVALMA.GT.NGMAXY) GOTO 50
  145. NA=IMMM(/1)
  146. IPRELL=IPREL
  147. DO 611 IA=1,NA
  148. C ... Si l'inconnue présente un mode d'ensemble ...
  149. IF(IMMM(IA).NE.0) THEN
  150. C ... On incrémente le compteur et
  151. IAA=IAA+1
  152. C ... on met le N° de la ligne dans ITRAA ...
  153. ITRAA(IAA)=IPRELL+IA-1
  154. ENDIF
  155. 611 CONTINUE
  156. C ... IFIB sera vu dans MONDE2 par le COMMON/CMOND2/ ...
  157. IFIB=IVPO(/1)
  158. segact lign
  159. CALL MONDE2(ITTRV(1),IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  160. > NA,IPREL,MULRE,INC,IVS,LLOM,IFIB,dnorma)
  161. 610 CONTINUE
  162. C ... Si on n'a pas eu de pb de mémoire on passe par là ...
  163. SEGSUP BID
  164. C ... ILMAX vaut le dernier noeud qui a pu être traité ...
  165. ILMAX=NNOE
  166. C ... On va à la division par la diagonale ...
  167. GOTO 55
  168.  
  169. 50 CONTINUE
  170. C ... Si on est là, c'est à cause des pb de mémoire ...
  171. SEGSUP BID
  172. C ... IVS est le N° du LIGN qui n'a pas pu être traité ...
  173. C ... ILMAX = N° du dernier traité ...
  174. ILMAX=IVS-1
  175. C ... IILMAX = N° du premier à traiter ...
  176. IILMAX=IVS
  177. DO 54 IVS=IILMAX,NNOE
  178. 58 CONTINUE
  179. LIGN=ILIGN(IVS)
  180. C ... Même si tout à l'heure SEGACT n'a pas marché, maintenant on a supprimé BID ...
  181. SEGACT /ERR=56/ LIGN
  182. GOTO 57
  183. 56 CONTINUE
  184. * EN CAS DE PROBLEME ON FAIT UN PEU DE PLACE
  185. C ... Si on a toujours des pb de mémoire, on désactive le LIGN N° ILMAX,
  186. LIGN=ILIGN(ILMAX)
  187. SEGDES LIGN*(NOMOD,MRU)
  188. C ... puis ce dernier est décrémenté ...
  189. ILMAX=ILMAX-1
  190. IF (ILMAX.EQ.1) CALL ERREUR(5)
  191. C ... et on recommence ...
  192. GOTO 58
  193. C ... On est là si tout s'est bien passé avec SEGACT ...
  194. 57 CONTINUE
  195. NA=IMMM(/1)
  196. IPRELL=IPREL
  197. DO 612 IA=1,NA
  198. IF(IMMM(IA).EQ.0) GOTO 612
  199. IAA=IAA+1
  200. ITRAA(IAA)=IPRELL+IA-1
  201. 612 CONTINUE
  202. IFIB = IVPO(/1)
  203. CALL MONDE2(ITTRV(1),IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  204. > NA,IPREL,MULRE,INC,IVS,LLOM,IFIB,dnorma)
  205. SEGDES,LIGN*(NOMOD,MRU)
  206. 54 CONTINUE
  207. 55 CONTINUE
  208. C
  209. C ... À cet endroit la descente est finie. L'état des LIGN est le suivant :
  210. C ... Ceux de IDEP à ILMAX sont actifs, les autres (si ILMAX < NNOE) sont désactivés ...
  211. C
  212. C **** DIVISION PAR LE TERME DIAGONAL ****
  213. C
  214. dnorma=0.d0
  215. idenorm=0
  216. DO 120 K=0,MUINC,INC
  217. DO 12 I=1,INC
  218. J=I+K
  219. if (2*abs(diag(i)).gt.abs(diag(i))) goto 122
  220. idenorm=1
  221. 122 continue
  222. VECTBB(J)=VECTBB(J)/DIAG(I)
  223. dnorma=max(dnorma,abs(vectbb(j)))
  224. 12 CONTINUE
  225.  
  226. 120 CONTINUE
  227. dnorma=dnorma*xzprec*xzprec
  228. SEGDES,MDIAG
  229.  
  230. C
  231. C **** MONTEE ****
  232. C
  233. C ... Si la matrice est asymétrique ...
  234. IF(EXISTS) THEN
  235.  
  236. C ... On commence par désactiver les LIGN qui sont restés actifs ...
  237. DO 2000 I=IDEP,ILMAX
  238. LIGN=ILIGN(I)
  239. SEGDES,LIGN*(NOMOD,MRU)
  240. 2000 CONTINUE
  241.  
  242. C ... Puis on désactive MILIGN ...
  243. SEGDES,MILIGN
  244.  
  245. C ... On change de MILIGN ...
  246. MILIGN=IILIGS
  247. SEGACT,MILIGN
  248.  
  249. C ... Et on active des LIGN pointés par ce dernier ...
  250. ILMAX=IDEP-1
  251. DO 2100 I=IDEP,NNOE
  252. LIGN=ILIGN(I)
  253. SEGACT /ERR=2110/ LIGN
  254. ILMAX = I
  255. 2100 CONTINUE
  256. C ... Et on passe à la montée ...
  257. GOTO 3000
  258.  
  259. C ... Si on n'a même pas activé le premier, alors ADIOS ...
  260. 2110 IF(ILMAX.EQ.IDEP-1) CALL ERREUR(5)
  261.  
  262. 3000 CONTINUE
  263. ENDIF
  264.  
  265.  
  266. J=NNOE
  267. C ... Début de la pseudo boucle sur les LIGN qui sont restés désactivés ...
  268. 70 CONTINUE
  269. IF (J.EQ.ILMAX) GOTO 72
  270. LIGN=ILIGN(J)
  271. SEGACT /ERR=68/ LIGN
  272. GO TO 67
  273. C ... Si on a des pb avec activation, on désactive et on diminue ILMAX ...
  274. 68 CONTINUE
  275. LIGN = ILIGN( ILMAX)
  276. SEGDES LIGN*(NOMOD,MRU)
  277. ILMAX=ILMAX-1
  278. IF(ILMAX.EQ.1) CALL ERREUR (5)
  279. C ... puis on recommence la tentative ...
  280. GO TO 70
  281. C ... On a réussi à activer ...
  282. 67 CONTINUE
  283. NA=IMMM(/1)
  284. IPRELL=IPREL
  285. IFIB=IVPO(/1)
  286. CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  287. > NA,IPREL,MULRE,INC,IFIB,dnorma)
  288. SEGDES,LIGN*(NOMOD,MRU)
  289. J = J-1
  290. GO TO 70
  291. 72 CONTINUE
  292. CC FIN DE PSEUDO BOUCLE J = INC ,ILMAX+1,-1 - Vieux commentaire et Faux !!!
  293. CC FIN DE PSEUDO BOUCLE J = NNOE ,ILMAX+1,-1
  294.  
  295. C ... Dans cette boucle on parcourt les LIGN qui sont restés actifs ...
  296. DO 13 J=ILMAX,IDEP,-1
  297. LIGN=ILIGN(J)
  298. NA=IMMM(/1)
  299. IPRELL=IPREL
  300. IFIB=IVPO(/1)
  301. CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  302. > NA,IPREL,MULRE,INC,IFIB,dnorma)
  303. SEGDES,LIGN*(NOMOD,MRU)
  304. 13 CONTINUE
  305.  
  306. C ... On n'avait même pas essayé d'activer les IDEP-1 premiers LIGN ...
  307. DO 113 J=IDEP-1,1,-1
  308. LIGN=ILIGN(J)
  309. SEGACT LIGN
  310. NA=IMMM(/1)
  311. IPRELL=IPREL
  312. DO 1131 ILM=1,NA
  313. IF(IMMM(ILM).EQ.0) GO TO 1131
  314. IAA=IAA+1
  315. ITRAA(IAA)=IPRELL+ILM-1
  316. 1131 CONTINUE
  317. IFIB=IVPO(/1)
  318. CALL MONDE1(IPPVV(1),VECTBB(1),VAL(1),IVPO(1),
  319. > NA,IPREL,MULRE,INC,IFIB,dnorma)
  320. SEGDES,LIGN*(NOMOD,MRU)
  321. 113 CONTINUE
  322.  
  323. if (idenorm.eq.1) then
  324. call erreur(1049)
  325. do k=0,muinc,inc
  326. do i=1,inc
  327. vectbb(i+k)=0.d0
  328. enddo
  329. enddo
  330. endif
  331. C
  332. C **** ON VERIFIE QUE PAS DE MODE RIGIDE ACTIVE
  333. C
  334. C ... IAA = nombre de modes d'ensemble (selon les indications dans IMMM) ...
  335. NENS=IAA
  336. NBEN=0
  337. IF(NENS.EQ.0) GO TO 26
  338.  
  339. MINCPO=IINCPO
  340. MIMIK=IIMIK
  341. SEGACT,MINCPO,MIMIK
  342. MELEME=IGEOMA
  343. SEGACT,MELEME
  344. NNOE=INCPO(/2)
  345. IINC1=INCPO(/1)
  346.  
  347. C ... Boucle sur les seconds membres ...
  348. DO 300 KI=0,MUINC,INC
  349. C ... XMA et XMAL seront le maximum des val. abs. des termes
  350. C du résultat partiel (avant la multiplication par MDNOR)
  351. C N° KI+1 multiplié par 1.e-10 ...
  352. XMA=xpetit/xzprec
  353. XMAL=xpetit/xzprec
  354. DO 30 KK=1,INC
  355. I=KK+KI
  356. if (ittr(kk).eq.0) then
  357. XMA=MAX(XMA,ABS(VECTBB(I)))
  358. else
  359. XMAL=MAX(XMAL,ABS(VECTBB(I)))
  360. endif
  361. 30 CONTINUE
  362. XMA=XMA*1.d-10
  363. XMAL=XMAL*1.d-10
  364. xmal=max(xma*1d-2,xmal)
  365. * write (6,*) ' mondes xma xmal ',xma,xmal
  366.  
  367. C ... Boucle sur les modes d'ensemble ...
  368. * write (6,*) ' nombre de modes d ensemble',nens
  369. iwrite = 0
  370. DO 21 IA=1,NENS
  371. I1=ITRAA(IA)
  372. J=IPNO(I1)
  373. DO 22 K=1,IINC1
  374. IF(INCPO(K,J).EQ.I1) GO TO 23
  375. 22 CONTINUE
  376. CALL ERREUR(5)
  377. RETURN
  378.  
  379. 23 CONTINUE
  380. C ... Si ce n'est pas un multiplicateur, le déplacement doit être
  381. C < à XMA, sinon ERREUR ...
  382. * write (6,*) ' mondes i1 ittr vect xma ',
  383. * > i1,ittr(i1),vectbb(i1+ki),xma
  384. IF(ITTR(I1).EQ.0.AND.ABS(VECTBB(I1+KI)).LE.XMA ) then
  385. vectbb(i1+ki)=0.d0
  386. GO TO 20
  387. endif
  388. C ... Si c'est un multiplicateur, le multiplicateur doit être
  389. C < à XMAL sinon ERREUR ...
  390. if (ittr(i1).ne.0) then
  391. i2=ittr(i1)
  392. if (abs(vectbb(i1+ki)+vectbb(i2+ki)).le.xmal) then
  393. vectbb(i1+ki)=0.d0
  394. vectbb(i2+ki)=0.d0
  395. goto 20
  396. endif
  397. endif
  398. * write (6,*) ' ittr vect ',ittr(i1),vectbb(i1+ki),
  399. * > vectbb(i2+ki),xmal
  400.  
  401.  
  402. C ... Si option NOEN on ne fait pas d'erreur, on accumule les pts dans un meleme
  403. C ... Si option NOEN on accumule les pts dans un meleme
  404. IF(NOEN.EQ.0) THEN
  405. nbelem=ipt8.num(/2)+1
  406. segadj ipt8
  407. ipt8.num(1,nbelem)=num(1,j)
  408. ENDIF
  409. NBEN=NBEN+1
  410. MOTERR(1:4)=IMIK(K)
  411. INTERR(1)=NUM(1,J)
  412. IF(NOEN.EQ.0) THEN
  413. GO TO 21
  414. ENDIF
  415. ITY=149
  416. CALL ERREUR(ITY)
  417. RETURN
  418.  
  419. 20 CONTINUE
  420. C ... Messages d'information ...
  421. JJK=NUM(1,J)
  422. KIK=KI/INC +1
  423.  
  424.  
  425. * on n'ecrit qu'une seule fois le message indetermination
  426. IF(ITTR(I1).EQ.0 .AND. MULRE.EQ.1.and.iwrite.eq.0) then
  427. if (imik(k).ne.'LX ')
  428. > WRITE(IOIMP,41) JJK,IMIK(K)
  429. iwrite=iwrite+1
  430. endif
  431.  
  432. IF (IIMPI.NE.0.AND. ITTR(I1).NE.0.AND.MULRE.EQ.1) then
  433. WRITE(IOIMP,42) JJK,IMIK(K)
  434. endif
  435.  
  436. IF(ITTR(I1).eq.0 .AND. MULRE.NE.1.and.iwrite.eq.0) then
  437. WRITE(IOIMP,410)KIK,JJK,IMIK(K)
  438. iwrite=iwrite+1
  439. endif
  440. IF (IIMPI.NE.0 .AND. ITTR(I1).NE.0.AND.MULRE.NE.1)
  441. & WRITE(IOIMP,420)KIK,JJK,IMIK(K)
  442.  
  443. 21 CONTINUE
  444.  
  445. 300 CONTINUE
  446.  
  447. 41 FORMAT(' INDETERMINATION DETECTEE AU NOEUD ',I6,' INCONNUE ',
  448. * A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ',
  449. * 'LA SUSDITE dans mondes')
  450. 42 FORMAT(' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ',
  451. * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ',
  452. * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE dans mondes')
  453. 410 FORMAT(' VECTEUR NUMERO',I3,' INDETERMINATION DETECTE AU NOEUD '
  454. *,I6,' INCONNUE ',
  455. * A4,/,' INDETERMINATION LEVEE PAR LA MISE A ZERO DE ',
  456. * 'LA SUSDITE dans mondes')
  457. 420 FORMAT(' VECTEUR NUMERO ',I3,/,
  458. * ' INDETERMINATION ENTRE CONDITIONS IMPOSEES DETECTEE ',
  459. * 'AU NOEUD ',I6,' INCONNUE ',A4,/,' INDETERMINATION LEVEE ',
  460. * 'PAR LA SUPPRESSION DE LA CONDITION REDONDANTE dans mondes')
  461.  
  462. SEGDES,MELEME
  463. SEGDES,MINCPO
  464. SEGDES,MIMIK
  465. 26 CONTINUE
  466. if (noen.eq.0) then
  467. * test si nan ou inf dans le resultat
  468. inan = 0
  469. DO 500 KI=0,MUINC,INC
  470. DO 500 KK=1,INC
  471. i = kk + ki
  472. if (abs(vectbb(i)).lt.xgrand) goto 500
  473. inan = inan + 1
  474. vectbb(i)=xgrand
  475. 500 continue
  476. if (inan.ne.0) nben = -inan
  477. endif
  478. * indiquer si besoin le nombre de modes d'ensemble excites
  479. * et le maillage des noeuds freinés
  480. IF (NOEN.EQ.0) then
  481. * write (6,*) ' mondes nben ',nben
  482. CALL ECRENT(NBEN)
  483. call ecrobj('MAILLAGE',ipt8)
  484. segdes ipt8
  485. endif
  486. MDNOR=IDNORM
  487. SEGACT,MDNOR
  488. DO 350 K=0,MUINC,INC
  489. DO 35 I = 1, INC
  490. VECTBB(I+K)=VECTBB(I+K)*DNOR(I)
  491. 35 CONTINUE
  492. * on force l'egalite des multiplicateurs de lagrange
  493. DO 36 I = 1, INC
  494. if (ITTR(I).ne.0) then
  495. * write (6,*) ' dans mondes ',i,ittr(i)
  496. if (vectbb(i+k).eq.0.d0.or.vectbb(ittr(i)+k).eq.0.d0) then
  497. * write (6,*) ' mondes vectbbs ',vectbb(i+k),vectbb(ittr(i)+k)
  498. vectbb(i+k)=0.d0
  499. vectbb(ittr(i)+k)=0.d0
  500. goto 36
  501. endif
  502. vectbb(i+k)=(vectbb(i+k)+vectbb(ittr(i)+k))/2
  503. vectbb(ittr(i)+k)=vectbb(i+k)
  504. endif
  505. 36 CONTINUE
  506. 350 CONTINUE
  507. SEGDES,MDNOR
  508. C ... On ne désactive MDNO1 que dans le cas où il est
  509. C différent de MDNOR ...
  510. IF(IDNORD.GT.0) THEN
  511. SEGDES,MDNO1
  512. ENDIF
  513. SEGDES MMATRI
  514. SEGDES,MILIGN
  515. SEGDES,MVECTD
  516. SEGSUP ITTRV
  517. IF(NE1.NE.0) SEGSUP,ITRAA
  518.  
  519. IF(IIMPI.EQ.3) WRITE(IOIMP,1001) MVECTD
  520. 1001 FORMAT(' SUBROUTINE MONDES : POINTEUR DU VECTEUR EN SORTIE=',I5)
  521.  
  522. CALL OOOMRU(0)
  523. RETURN
  524. END
  525.  
  526.  
  527.  

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