Télécharger ldmt2.eso

Retour à la liste

Numérotation des lignes :

ldmt2
  1. C LDMT2 SOURCE PV 21/04/26 21:15:10 10866
  2. SUBROUTINE LDMT2(ITRAV1,ITOPO1,INUIN1,IMINI1,MMMTRI,IPO1,INCTR1
  3. & ,INCTR2,IITOP1,TRSUP,MENAGE,LDIAG,IITOPB,
  4. & ITOPOB,IPODD,njtot,INORMU)
  5.  
  6. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  7. C
  8. C **** SUBROUTINE POUR FAIRE L'ASSEMBLAGE DE MATRICES
  9. C
  10. C EN ENTREE:
  11. C **** ITRAV1 : POINTEUR OBJET MRIGIDITE
  12. C **** ITOPO1 : POINTEUR SEGMENT DE TRAVAIL ITOPO ( VOIR ASSEM1)
  13. C **** IITOP1 : POINTEUR SEGMENT DE TRAVAIL IITOP ( VOIR ASSEM1)
  14. C **** INUIN1 : POINTEUR SEGMENT DE TRAVAIL INUINV(VOIR ASSEM1)
  15. C **** IMINI1 : POINTEUR SEGMENT DE TRAVAIL IMINI (VOIR ASSEM1)
  16. C **** IPO1 : POINTEUR SEGMENT DE TRAVAIL IPOS (VOIR ASSEM1)
  17. C **** MMMTRI : POINTEUR OBJET MATRICE TRIANGULARISEE (NON MODIFIE)
  18. C (VOIR SMMATRI)
  19. C **** INORMU : =0 si on ne veut pas normaliser les multiplicateurs
  20. C de Lagrange, <>0 sinon
  21. C
  22. C Appelée par : LDMT1
  23. C
  24. C Auteur : Michel BULIK
  25. C
  26. C Date : Printemps '95
  27. C
  28. C Langage : ESOPE + FORTRAN77
  29. C
  30. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  31. IMPLICIT INTEGER(I-N)
  32. IMPLICIT REAL*8 (A-H,O-Z)
  33. -INC CCREEL
  34.  
  35. -INC PPARAM
  36. -INC CCOPTIO
  37. -INC SMELEME
  38. -INC SMRIGID
  39. -INC SMMATRI
  40.  
  41. SEGMENT,INUINV(NNGLOB)
  42. SEGMENT,ITOPO(IENNO)
  43. SEGMENT,ITOPOB(IENNO)
  44. SEGMENT,IITOP(NNOE+1)
  45. SEGMENT,IITOPB(NNOE+1)
  46. SEGMENT,IMINI(INC)
  47. SEGMENT,IPOS(NNOE1)
  48. SEGMENT,IPOD(NNOE1)
  49. SEGMENT,INCTRR(NIRI)
  50. SEGMENT,INCTRD(NIRI)
  51. SEGMENT,INCTRA(NLIGRE)
  52. SEGMENT,INCTRB(NLIGRE)
  53. SEGMENT,IPV(NNOE)
  54. SEGMENT,VMAX(INC)
  55. C
  56. C **** CES TABLEAUX SERVENT AU REPERAGE DE LA MATRICE POUR L'ASSEMBLAG
  57. C **** IL SERONT TOUS SUPPRIMES EN FIN D'ASSEMBLAGE.
  58. C
  59. SEGMENT,IVAL(NNN)
  60. SEGMENT,ITRA(NNN,2)
  61. SEGMENT TRATRA
  62. REAL*8 XTRA(INCRED,INCDIF)
  63. INTEGER LTRA(INC,INCDIF)
  64. INTEGER NTRA(INCRED,INCDIF)
  65. INTEGER MTRA(INCDIF)
  66. ENDSEGMENT
  67.  
  68. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  69. C
  70. C **** IVAL(I)=J : LA I EME LIGNE D'UNE PETITE MATRICE S'ASSEMBLE
  71. C DANS LA J EME DE LA GRANDE.
  72. C **** ITRAV(I,1)=J : LA IEME INCONNUE DU NOEUD EN COURS D'ASSEMBLAGE
  73. C ET QUI SE TROUVE DANS LA PETITE MATRICE SE TROUVE
  74. C EN J EME POSITION DE LA PETITE MATRICE.
  75. C **** ITRAV(I,2) : LA IEME INCONNUE DU NOEUD EN COURS D'ASSEMBLAGE
  76. C PRESENT DANS LA PETITE MATRICE EST EN JEME
  77. C POSITION DANS LA GRANDE
  78. C
  79. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  80.  
  81. SEGMENT,RA(N1,N1)*D
  82. SEGMENT JNOMUL
  83. LOGICAL INOMUL(NNR)
  84. ENDSEGMENT
  85. REAL*8 DMAX,COER,DDDD,DMAXY
  86. LOGICAL NOMUL
  87. LOGICAL TRSUP,TEST,MENAGE,LDIAG
  88.  
  89.  
  90. SAVE IPASG
  91. DATA IPASG/0/
  92.  
  93.  
  94. ** write(6,*) ' ldmt2 inormu ',inormu
  95. C
  96. C en cas de normalisation des variables.
  97. C
  98. C ... On a vérifié dans RESOU que dans le cas non-symétrique soit NORINC=NORIND=0,
  99. C soit les deux sont non-nuls, mais on se méfie quand même ...
  100. inwuit=0
  101. C write (6,*) ' NORINC NORIND LDIAG',norinc,norind,ldiag
  102. IF(NORINC.NE.0.AND.NORIND.NE.0 .AND. LDIAG) then
  103. CALL ASSEM0(ITRAV1,1,inwuit)
  104. endif
  105. C
  106. C PV ON ACTIVE UNE FOIS POUR TOUTES LES MELEME DESCR... DE LA RIGIDITE
  107. C ON EN PROFITE POUR CREER INOMUL
  108. C
  109. C **** RECHERCHE DE LA DIMENSION MAX DE IVAL,ET SEGINI DE IVAL ET ITRA
  110. C
  111. INCTRR=INCTR1
  112. SEGACT,INCTRR
  113. C ... dans le cas asymétrique on se doit de faire la même chose avec les vd ...
  114. INCTRD=INCTR2
  115. SEGACT,INCTRD
  116.  
  117. MRIGID=ITRAV1
  118. SEGACT,MRIGID
  119. NNR=IRIGEL(/2)
  120. C ... NNR = équivalent à NRIGEL (nb. de rigidités élémentaires) ...
  121.  
  122. NNN=0
  123. SEGINI JNOMUL
  124. DO 1 IRI=1,NNR
  125.  
  126. INCTRA=INCTRR(IRI)
  127. SEGACT INCTRA
  128. C ... dans le cas asymétrique on se doit de faire la même chose avec les vd ...
  129. INCTRB=INCTRD(IRI)
  130. SEGACT INCTRB
  131.  
  132. DESCR=IRIGEL(3,IRI)
  133. SEGACT, DESCR
  134.  
  135. IPT1=IRIGEL(1,IRI)
  136. SEGACT IPT1
  137.  
  138. xMATRI=IRIGEL(4,IRI)
  139. SEGACT xMATRI
  140.  
  141. if (.not.trsup) then
  142. C ... NA = nb de variables primales ...
  143. NA=LISINC(/2)
  144. else
  145. C ... NA = nb de variables duales ...
  146. NA=LISDUA(/2)
  147. endif
  148. C ... NNN = max des nb de variables primales ...
  149. NNN=MAX(NA,NNN)
  150.  
  151. C ... INOMUL (NOn MULtiplicateurs ?) dit s'il ne s'agit pas des CL ...
  152. INOMUL(IRI)=.TRUE.
  153. IF(IPT1.ITYPEL.EQ.49 )INOMUL(IRI)=.FALSE.
  154.  
  155. 1 CONTINUE
  156. C ... IVAL a la taille NNN ...
  157. SEGINI,IVAL
  158. C ... ITRA A LA TAILLE (NNN,2) ...
  159. SEGINI,ITRA
  160. C
  161. C **** ACTIVATION DES SEGMENTS DE TRAVAILS ET DE MMATRI
  162. C
  163. IF(.NOT.LDIAG) THEN
  164. ITOPO=ITOPO1
  165. IITOP=IITOP1
  166. IPOS=IPO1
  167. ELSE
  168. ITOPO=ITOPOB
  169. IITOP=IITOPB
  170. IPOS = IPODD
  171. ENDIF
  172. SEGACT,IITOP
  173. SEGACT,ITOPO
  174. SEGACT,IPOS
  175. C write(6,*) 'ipos', ( ipos(IU),IU=1,ipos(/1))
  176. C write(6,*) 'itopo', ( itopo(IU),IU=1,itopo(/1))
  177. C write(6,*) 'iitop', ( iitop(IU),IU=1,iitop(/1))
  178. INUINV=INUIN1
  179. SEGACT,INUINV
  180. SEGACT,IPOS
  181. C ... NNOE = nombre de noeuds impliqués ...
  182. NNOE=IPOS(/1)-1
  183. C ... INC = dimension de la matrice ...
  184. INC=IPOS(NNOE+1)
  185. MMATRI=MMMTRI
  186. SEGACT,MMATRI*MOD
  187. IF(LDIAG) THEN
  188. SEGINI,MDIAG
  189. IDIAG=MDIAG
  190. ELSE
  191. MDIAG=IDIAG
  192. SEGACT,MDIAG
  193. ENDIF
  194. SEGINI,MILIGN
  195. IF(TRSUP) THEN
  196. IILIGS=MILIGN
  197. ELSE
  198. IILIGN=MILIGN
  199. ENDIF
  200. MINCPO=IINCPO
  201. SEGACT,MINCPO
  202. C ... INCDIF = nb de variables primales différentes ...
  203. INCDIF=INCPO(/1)
  204. MIPO1=IDUAPO
  205. SEGACT,MIPO1
  206. casym(OK) ... on va prendre en compte aussi les vd pour déterminer
  207. c la taille de TRATRA ...
  208. INCDID=MIPO1.INCPO(/1)
  209. INCDIF=MAX(INCDIF,INCDID)
  210.  
  211. MIMIK=IIMIK
  212. SEGACT,MIMIK
  213.  
  214. C ... IPV est de taille NNOE ...
  215. SEGINI IPV
  216.  
  217. INCRED=0
  218. DO 80 INO =1,NNOE
  219. ICOMPT=0
  220. C ... MAXELE = nb d'éléments qui contiennent le noeud No INO ...
  221. MAXELE = (IITOP(INO+1)-IITOP(INO))/2
  222. C ... Boucle sur ces éléments ...
  223. DO 81 IELE=1,MAXELE
  224.  
  225. C ... IIU = pointeur dans IITOP ...
  226. IIU=IITOP(INO) + IELE + IELE -2
  227.  
  228. C ... IEL = numero de l'élément dans son maillage ...
  229. IEL=ITOPO(IIU)
  230.  
  231. C ... IRI = numero du maillage (dans IRIGEL) de cet élément ...
  232. IRI=ITOPO(IIU+1)
  233.  
  234. C ... On va chercher ce maillage dans IRIGEL ...
  235. MELEME=IRIGEL(1,IRI)
  236.  
  237. C ... Boucle sur les noeuds de cet élément ...
  238. DO 83 I=1,NUM(/1)
  239.  
  240. C ... IP = numéro local du I-ème noeud de l'élément ...
  241. IP=INUINV(NUM(I,IEL))
  242.  
  243. C ... La première condition = les "connections" sont symétriques ...
  244. IF (IP.GT.INO) GOTO 83
  245. C ... Le noeud IP a déjà vu ce noeud ...
  246. IF (IPV(IP).EQ.INO) GOTO 83
  247.  
  248. C ... Ces deux opérations se font si (IP <= INO et IPV(IP) != INO) ...
  249. C ... IPV = N° max des noeuds "connectés" au noeud IP ...
  250.  
  251. IPV(IP)=INO
  252. C ... ICOMPT = nombre de noeuds connectés (de numéro < ou =) au noeud INO ...
  253. ICOMPT=ICOMPT+1
  254. 83 CONTINUE
  255. 81 CONTINUE
  256. C ... INCRED = maximum des ICOMPT sur tous les noeuds impliqués ...
  257.  
  258. INCRED=MAX(INCRED,ICOMPT)
  259. cdbg write(*,*) 'LDMT2 : Noeud ',INO,', ICOMPT = ',ICOMPT,
  260. cdbg & ' -> INCRED = ',INCRED
  261. 80 CONTINUE
  262.  
  263. cdbg segprt,ipv
  264. SEGSUP IPV
  265.  
  266. INCRED=INCRED*INCDIF
  267. C ... La taille de TRATRA dépend de INC, INCRED et INCDIF ...
  268. SEGINI TRATRA
  269. segini vmax
  270. C
  271. C ... Coefficients de normalisation ...
  272. C -----------------------------
  273.  
  274. LLVNUL=0
  275. C ... On n'initialise plus IJMAX ici - ceci se fait un niveau au dessus ...
  276. C-OLD IJMAX=0
  277.  
  278. C ... La taille de MDNOR est INC (nombre d'équations) ...
  279. MDNOR=0
  280. IF(LDIAG) THEN
  281. SEGINI MDNOR
  282. IDNORM=MDNOR
  283. C ... On a vérifié dans RESOU que dans le cas non-symétrique soit NORINC=NORIND=0,
  284. C soit les deux sont non-nuls ...
  285. SEGINI,MDNO1
  286. IDNORD=MDNO1
  287. DO 20 I=1,INC
  288. DNOR(I)=1.D0
  289. MDNO1.DNOR(I)=1.D0
  290. 20 CONTINUE
  291. else
  292. mdnor=idnorm
  293. segact mdnor
  294. mdno1=idnord
  295. segact mdno1
  296. ENDIF
  297.  
  298. C ... Si une normalisation est imposée, les valeurs dans MDNOR et MDNO1 seront modifiées par ASNS3 ...
  299. midua=iidua
  300. IF(NORINC.NE.0 .AND. NORIND.NE.0 .AND. LDIAG)
  301. & CALL ASNS3(MDNOR,MIMIK,MINCPO,mdno1,midua,mipo1)
  302.  
  303. C **** BOUCLE *100* SUR LES NUMEROS DE NOEUDS QUE L'ON ASSEMBLE
  304. C -------------------------------------------------------------
  305.  
  306. DO 100 INO=1,NNOE
  307. cdbg write(*,*) 'LDMT2: ----------------------------------'
  308. cdbg write(*,*) 'LDMT2: Noeud N° ',INO
  309. DO 101 IIT=1,INCDIF
  310. MTRA(IIT)=0
  311. 101 CONTINUE
  312.  
  313. C ... Attention à l'égalité des IPOS et IPOD qui est testée dans ASNS1 ...
  314.  
  315. C ... IPRE = N° de la première ligne (ou colonne) concernée par le noeud INO ...
  316. IPRE=IPOS(INO )+1
  317. C ... IDER = N° de la dernière ligne (ou colonne) concernée par le noeud INO ...
  318. IDER=IPOS(INO+1)
  319. cdbg write(*,*) 'LDMT2: Lignes de ',IPRE,' jusqu''à ',IDER
  320. LLVVA=0
  321. C
  322. C **** BOUCLE *99* SUR LES ELEMENTS TOUCHANT LE NOEUD INO
  323. C POUR LES ELEMNTS MULTIPLICATEUR ON NE FAIT PAS
  324. C L'ASSEMBLAGE
  325. C
  326. C ... MAXELE = nb d'éléments qui contiennent le noeud No INO ...
  327. MAXELE= (IITOP(INO+1) -IITOP(INO))/2
  328. DO 99 IELE=1,MAXELE
  329. cdbg write(*,*) 'LDMT2: - - - - - - - - - - - -'
  330. cdbg write(*,*) 'LDMT2: Elément N° ',IELE
  331. IIU=IITOP(INO) + IELE + IELE - 2
  332. C ... IEL = numero de l'élément dans son maillage ...
  333. IEL=ITOPO(IIU)
  334. C ... IRI = numero du maillage (dans IRIGEL) de cet élément ...
  335. IRI=ITOPO(IIU+1)
  336. C ... MELEME = pointeur vers ce maillage ...
  337. MELEME=IRIGEL(1,IRI)
  338. C ... DESCR = pointeur vers le segment DESCR concerné ...
  339. DESCR=IRIGEL(3,IRI)
  340. C ... INCTRA contient les indices (dans IMIK et IHAR) des descriptions
  341. C des DDL correspondants ...
  342. casym(OK) ... ici il faudra choisir en fonction de TRSUP si prendre INCTRR ou INCTRD ...
  343. INCTRA=INCTRR(IRI)
  344. INCTRB=INCTRD(IRI)
  345. C ... xMATRI = tableau de pointeurs vers les matrices élémentaires ...
  346. xMATRI=IRIGEL(4,IRI)
  347. C ... COER = coefficient multiplicateur ...
  348. COER=COERIG(IRI)
  349. C ... NOn MULtiplicateur ? ...
  350. NOMUL=INOMUL(IRI)
  351. C
  352. C **** NOMUL =.FALSE. IL EXISTE UN MULTIPLICATEUR
  353. C **** INITIALISATION DE IVAL. IVAL(I)=J VEUT DIRE QUE
  354. C **** LA I EME LIGNE DE LA PETITE MATRICE S'ASSEMBLE DANS
  355. C **** LA J EME DE LA GRANDE MATRICE.
  356. C
  357. C ... Maintenant, dans la boucle 96 (nouvelle par rapport à ASSEM2) on
  358. C va parcourir les lignes (colonnes) pour trouver celles, qui
  359. C s'appuyent sur le noeud INO, cas inférieur (supérieur) ...
  360.  
  361. NA=0
  362. IF(TRSUP) THEN
  363. NIN=LISINC(/2)
  364. ELSE
  365. NIN=LISDUA(/2)
  366. ENDIF
  367. DO 96 ICO=1,NIN
  368. IF(TRSUP) THEN
  369. IJA=INUINV(NUM(NOELEP(ICO),IEL))
  370. IJB=INCTRA(ICO)
  371. ELSE
  372. IJA=INUINV(NUM(NOELED(ICO),IEL))
  373. IJB=INCTRB(ICO)
  374. ENDIF
  375. C ... Si on a trouvé le bon noeud ...
  376. IF(IJA.EQ.INO) THEN
  377. C ... On incrémente NA (le nombre de DDL connectés au noeud INO ?) ...
  378. NA=NA+1
  379. cdbg write(*,*) 'LDMT2: IJA == INO => NA devient ',NA
  380. ITRA(NA,1)=ICO
  381. IF(TRSUP) THEN
  382. ITRA(NA,2)=INCPO(IJB,IJA)
  383. ELSE
  384. ITRA(NA,2)=MIPO1.INCPO(IJB,IJA)
  385. ENDIF
  386. ENDIF
  387. 96 CONTINUE
  388.  
  389. C ... Dans la boucle 98 on va remplir le tableau IVAL avec les numéros
  390. C de colonnes (lignes) pour tous les DDLs primaux (duaux) de la matrice,
  391. C ceci concerne le triangle inférieur (supérieur) ...
  392. casym(OK) ... la boucle devra se faire soit sur les duales soit sur les primales ...
  393.  
  394. IF(TRSUP) THEN
  395. NIN=LISDUA(/2)
  396. ELSE
  397. NIN=LISINC(/2)
  398. ENDIF
  399. C ... Boucle sur les variables primales (duales) de l'élément ...
  400. DO 98 ICO=1,NIN
  401. IF(TRSUP) THEN
  402. IJA=INUINV(NUM(NOELED(ICO),IEL))
  403. IJB=INCTRB(ICO)
  404. IVAL(ICO)=MIPO1.INCPO(IJB,IJA)
  405. cdbg write(*,*) 'LDMT2: vd locale N°',ICO,'-> vd globale N°',IJB
  406. cdbg write(*,*) 'LDMT2: => N° d''équation = ',IVAL(ICO)
  407. ELSE
  408. C ... IJA = numéro local du noeud-support du DDL ...
  409. casym(OK) ... ici il faudra choisir en fonction de TRSUP si prendre NOELEP ou NOELED ...
  410. IJA=INUINV(NUM(NOELEP(ICO),IEL))
  411. C ... IJB = indice du DDL ...
  412. casym(OK) ... INCTRA ou INCTRB ??? ...
  413. IJB=INCTRA(ICO)
  414. C ... On met dans IVAL(ICO) le N° de la colonne correspondant au noeud IJA
  415. C et DDL No IJB ...
  416. casym(OK) ... ici il faudra choisir en fonction de TRSUP si prendre MINCPO ou MIPO1 (?) ...
  417. IVAL(ICO)=INCPO(IJB,IJA)
  418. cdbg write(*,*) 'LDMT2: vp locale N°',ICO,'-> vp globale N°',IJB
  419. cdbg write(*,*) 'LDMT2: => N° de la colonne = ',IVAL(ICO)
  420. ENDIF
  421. 98 CONTINUE
  422.  
  423. C XMATRI=IMATTT(IEL)
  424. C SEGACT,XMATRI
  425. cdbg segprt,xmatri
  426. C
  427. C **** BOUCLE *95* SUR LES INCONNUES DE LA PETITE MATRICE
  428. C
  429. DO 95 INCC=1,NA
  430. C ... INCO = le N° de la ligne (cas inf.) ou colonne (cas sup.)
  431. C du DDL dual (primal) N° INCC du noeud INO ...
  432. casym(OK) ... ici il faudra choisir en fonction de TRSUP quoi prendre ...
  433. INCO=ITRA(INCC,2)
  434.  
  435. C ... Boucle sur les DDL's du noeud INO ...
  436. DO 90 IK=1,NIN
  437. C ... IO = N° de la colonne (la ligne) correspondant au DDL N° IK de l'élément
  438. casym(OK) ... ici il faudra choisir en fonction de TRSUP quoi prendre ...
  439. IO=IVAL(IK)
  440.  
  441. C ... Si IO > INCO on ne fait plus rien dans cette boucle ...
  442. C ... Car on ne parcourt que la moitié de la matrice !!! ...
  443. IF(IO.GT.INCO) GO TO 90
  444.  
  445. C ... ILOC = N° de ligne relatif / IPRE du noeud INO ...
  446. ILOC=INCO-IPRE+1
  447.  
  448. C ... JJ = N° (dans l'élément) de la ligne (de la colonne) correspondant au
  449. C DDL N° INCC du noeud INO ...
  450. casym(OK) ... attention à ce qui se passe ici ... !!!
  451. JJ=ITRA(INCC,1)
  452.  
  453. C ... Ci-dessous on stockera les quelques lignes de la rigidité concernant le
  454. C noeud INO dans le segment TRATRA. Les coefficients vont dans XTRA, le deuxième
  455. C indice (ILOC) est le numéro relatif de la ligne, le premier (IMTT) n'a rien à voir avec le
  456. C numéro de la colonne. Celui-ci est stocké dans NTRA(IMTT,ILOC). L'information inverse
  457. C (IMTT en fonction de IO) se trouve dans LTRA(IO,ILOC). On met dans MTRA le nombre
  458. C de termes différents assemblés par ligne ...
  459. C ... ILTT doit contenir l'indice IMMTT en fonction du N° de la colonne; s'il est
  460. C nul, le terme n'a pas encore été assemblé, il faudra donc initialiser ...
  461. ILTT= LTRA(IO,ILOC)
  462. IF(ILTT.EQ.0) THEN
  463. C ... LLVVA est initialisé à 0 au début de la boucle 100 et sert à compter les
  464. C termes assemblés ...
  465. LLVVA=LLVVA+1
  466. C ... Les MTRA sont mis à 0 au début de la boucle 100 ...
  467. C ... Ici on incrémente MTRA(ILOC) et on met la nouvelle valeur dans IMMTT,
  468. C qui sera le premier indice dans différents tableaux de TRATRA ...
  469. IMMTT=MTRA(ILOC)+1
  470. MTRA(ILOC)=IMMTT
  471. C ... puis on initialise le nouveau terme ...
  472. XTRA(IMMTT,ILOC)=0.D0
  473. C ... son N° de colonne (cas inf.) ou ligne (cas supérieur) ...
  474. NTRA(IMMTT,ILOC)=IO
  475. C ... et l'information inverse ...
  476. LTRA(IO,ILOC)=IMMTT
  477. ILTT=IMMTT
  478. ENDIF
  479. C ... Si ce n'est pas une condition aux limites, on assemble ...
  480. IF(NOMUL) THEN
  481. C ... Attention ! Dans XTRA 1er indice est lié à la colonne, deuxième est le numéro
  482. C de la ligne, tandis que dans RE c'est l'inverse ...
  483. IF(TRSUP) THEN
  484. C write(6,*) ' iltt iloc ' , iltt,iloc
  485. C write(6,*) 'ik jj iel re coer ',ik,jj,iel,re(ik,jj,iel),coer
  486. XTRA(ILTT,ILOC)=XTRA(ILTT,ILOC)+RE(IK,JJ,iel)*COER
  487. cdbg write(*,*) 'LDMT2: On rajoute le terme L',IK,', C',JJ,' soit',
  488. cdbg & RE(IK,JJ),' à la col',IPRE+ILOC-1,' et lig',
  489. cdbg & NTRA(ILTT,ILOC)
  490. ELSE
  491. C write(6,*) ' iltt iloc ' , iltt,iloc
  492. C write(6,*) 'ik jj iel re coer ',ik,jj,iel,re(jj,ik,iel),coer
  493. XTRA(ILTT,ILOC)=XTRA(ILTT,ILOC)+RE(JJ,IK,iel)*COER
  494. cdbg write(*,*) 'LDMT2: On rajoute le terme L',JJ,', C',IK,' soit',
  495. cdbg & RE(JJ,IK),' à la lig',IPRE+ILOC-1,' et col',
  496. cdbg & NTRA(ILTT,ILOC)
  497. ENDIF
  498.  
  499. ENDIF
  500. 90 CONTINUE
  501. 95 CONTINUE
  502. C SEGDES,XMATRI
  503. 99 CONTINUE
  504. C
  505. C *** COMPACTAGE DES LIGNES, EN MEME TEMPS CALCUL DE IJMAX QUI SERA
  506. C *** LA DIMENSION MAX D'UN SEGMENT LIGN.
  507. C *** LE SEGMENT ASSOCIE A UNE LIGNE (SEGMENT LLIGN) EST DE LA FORME :
  508. C *** IMMMM(NA) PERMET DE SAVOIR SI UN MOUVEMENT D'ENSEMBLE SUR LA
  509. C *** LIGNE EXISTE. IPPO(NA+1) DONNE LA POSITION DANS XXVA LA 1ERE
  510. C *** VALEUR DE LA LIGNE .XXVA VALEUR DE LA MATRICE.
  511. C *** LINC(I) DONNE LE NUMERO DE LA COLONNE DU IEME ELEM DE XXVA
  512. C
  513. C ... NA = nb de lignes concernant le noeud INO ...
  514. NA = IDER-IPRE+1
  515. C ... LLVNUL = somme cumulée des LLVVA ...
  516. LLVNUL=LLVNUL+LLVVA
  517. SEGINI,LLIGN
  518. ILIGN(INO)=LLIGN
  519. NBA=0
  520. C write(6,*) ' ider ipre llvnul ino ',ider,ipre,llvnul,ino
  521. C ... Boucle sur les DDL du noeud INO ...
  522. DO 120 JPA=1,NA
  523. C ... IIIN = N° de ligne du DDL N° JPA ...
  524. IIIN=IPRE+JPA -1
  525. C ... que l'on met au bon endroit dans IMMMM (faisant partie de LLIGN) ...
  526. IMMMM(JPA)=IIIN
  527. C ... IPPO(I)+1 = adresse du début (dans XXVA et LINC) des données relatives au DDL N° I ...
  528. IPPO(JPA)=NBA
  529. C ... Boucle sur les termes dans la ligne ...
  530. DO 121 IPAK = 1,MTRA(JPA)
  531. C ... IUNPAK = N° de la colonne du terme N° IPAK ...
  532. IUNPAK=NTRA(IPAK,JPA)
  533. C ... Pour les termes mis dans LLIGN on efface l'information sur la position dans XTRA ...
  534. LTRA(IUNPAK,JPA)=0
  535. C ... Compteur ++ ...
  536. NBA=NBA+1
  537. C ... Le N° de la colonne va dans LINC ...
  538. LINC(NBA)=IUNPAK
  539. C write(6,*) ' nba ipak jpa, ', nba ,ipak,jpa
  540. C ... Transfert du XTRA (segment TRATRA) vers XXVA (segment LLIGN) ...
  541. XXVA(NBA)=XTRA(IPAK,JPA)
  542. C ... Les termes diagonaux vont dans DIAG (segment MDIAG) ...
  543. IF(LDIAG) THEN
  544. C write(6,*) ' iiin nba xxva(nba)', iiin,nba,xxva(nba)
  545. IF(IIIN.EQ.IUNPAK) DIAG(IIIN)=XXVA(NBA)
  546. ENDIF
  547. vmax(iiin)=max(abs(xxva(nba)),vmax(iiin))
  548. 121 CONTINUE
  549. 120 CONTINUE
  550. C ... Le pointeur vers la fin de la dernière ligne ...
  551. IPPO(NA+1)= NBA
  552. NJMAX= 0
  553. C recherche du mini globale sur toutes les inconnues
  554. LPA=IPRE
  555. C ... Boucle sur les lignes relatives au noeud INO ...
  556. DO 126 JPA=IPRE,IDER
  557. C ... On met le N° du noeud dans IPNO (segment MILIGN) ...
  558. IPNO(JPA)=INO
  559. C ... IPDE et IPDF : début et fin des données relatives au DDL N° JPA dans XXVA et LINC ...
  560. IPDE=IPPO(JPA-IPRE+1)+1
  561. IPDF=IPPO(JPA-IPRE+2)
  562. C ... LPA = N° mini de la colonne avec des termes non nuls dans la ligne N° JPA ...
  563. DO 155 JHT=IPDE,IPDF
  564. LPA=MIN(LPA,LINC(JHT))
  565. 155 CONTINUE
  566. 126 CONTINUE
  567. DO 127 JPA=IPRE,IDER
  568. C ... On le met dans LDEB (segment LLIGN) ...
  569. LDEB(JPA-IPRE+1)=LPA
  570. C ... NNA = longueur de la "skyline" ...
  571. NNA= JPA- LPA +1
  572. C ... NJMAX = profil cumulé sur un noeud ...
  573. NJMAX=NJMAX+NNA
  574. 127 CONTINUE
  575. C ... NJTOT = profil total ...
  576. NJTOT=NJTOT+NJMAX
  577. C ... IJMAX = (profil / noeud) maxi ...
  578. IF(IJMAX.LT.NJMAX) IJMAX=NJMAX
  579. SEGDES,LLIGN
  580. 100 CONTINUE
  581. C ... Fin de la boucle sur les noeuds ...
  582. SEGSUP TRATRA
  583. C
  584. C **** ON REPREND TOUTE LES MATRICES CONTENANT LES MULTIPLICATEURS
  585. C **** POUR MULTIPLIER TOUS LEURS TERMES PAR UNE NORME ATTACHEE
  586. C **** A CHAQUE MULTIPLICATEUR.PUIS ON LES ASSEMBLE.
  587. C
  588. * d'abord etablir une norme generale pour le cas ou on n'arrive pas
  589. * a calculer la norme particuliere
  590. DMAXGE=xpetit
  591. DO 378 I=1,INC
  592. ** write (6,*) ' assem2 diag vmav ',diag(i),vmax(i)
  593. ** vmax(i)=max(vmax(i),abs(diag(i)))
  594. DMAXGE=MAX(DMAXGE,abs(vmax(i)))
  595. 378 CONTINUE
  596. if (iimpi.ne.0 )
  597. > write (6,*) ' nb inconnues facteur multiplicatif general ',
  598. > INC,DMAXGE
  599. if (dmaxge.lt.xpetit/xzprec*10) dmaxge=1.d0
  600. IENMU=0
  601. C ... ATTENTION ! Ici commence une boucle cachée (exécutée avec un GOTO 375) ...
  602. 375 IENMU1 = IENMU
  603. IENMU=0
  604. C ... Boucle sur les rigidités qui calcule le nombre de matrices de blocages ...
  605. DO 376 I=1,NNR
  606. IF(.NOT.INOMUL(I)) IENMU=IENMU+1
  607. 376 CONTINUE
  608.  
  609. C ... S'il n'y en a pas on saute cette partie du code ...
  610. IF( IENMU.EQ.0) GO TO 3750
  611.  
  612. C ... MIMIK contient les noms des variables primales ...
  613. MIMIK=IIMIK
  614. SEGACT,MIMIK
  615. C ... Boucle sur les rigidités ...
  616. DO 11 I=1,NNR
  617. C ... Si celle si n'est pas une matrice de blocage on passe à la suivante ...
  618. IF(INOMUL(I)) GO TO 11
  619.  
  620. DESCR=IRIGEL(3,I)
  621. C ... N3 = nb de variables primales ...
  622. N3=LISINC(/2)
  623. COER=COERIG(I)
  624. MELEME=IRIGEL(1,I)
  625. INCTRA=INCTRR(I)
  626. xMATRI=IRIGEL(4,I)
  627. C ... N2 = nombre d'éléments dans le support géométrique ...
  628. N2=NUM(/2)
  629. C ... À quoi correspond ce cas ? (pas de matrices élémentaires) ...
  630. IF (re(/3).EQ.0) THEN
  631. INOMUL(I)=.TRUE.
  632. GOTO 11
  633. ENDIF
  634. C XMATRI=IMATTT(1)
  635. C SEGACT,XMATRI
  636. C ... N1 = nombre de variables duales ...
  637. C ... Pourquoi va-t-on chercher N3 dans DESCR et N1 dans RE ? ...
  638. N1=RE(/1)
  639. SEGINI,RA
  640. C ... Boucle sur les éléments ...
  641. DO 14 IEL=1,N2
  642.  
  643. C ... Boucle sur les inconnues ...
  644. DO 15 ICO=1,N3
  645. C ... IJA = N° local du noeud-support de l'inconnue N° ICO ...
  646. IJA=INUINV(NUM(NOELEP(ICO),IEL))
  647. C ... IJB = N° de l'inconnue ...
  648. IJB=INCTRA(ICO)
  649. C ... IVAL contient la correspondance entre le N° local du DDL
  650. C et le N° d'équation (de ligne) correspondant ...
  651. IVAL(ICO)=INCPO(IJB,IJA)
  652. 15 CONTINUE
  653. C JCARDO => pour les fous qui veulent se passer de la
  654. C normalisation des mult. de Lagrange :)
  655. IF (INORMU.EQ.0) THEN
  656. DMAX=1.D0
  657. DMAXY=DMAX
  658. ELSE
  659. DDDD= xpetit
  660. C ... DMAX = max des valeurs absolues des termes diagonaux correspondants
  661. C au DDL de l'élément N° IEL ...
  662. DMAX=xpetit
  663. IF((.NOT.LDIAG).AND.(IDIAG.EQ.0)) THEN
  664. C write(*,*) 'Erreur interne dans LDMT2 : MDIAG pas encore',
  665. C & ' initialisé !!!'
  666. RETURN
  667. ENDIF
  668. DO 19 ICO=1,N3
  669. DMAX=MAX(DMAX,vmax(IVAL(ICO)),abs(diag(ival(ico))))
  670. 19 CONTINUE
  671.  
  672. C ... AUX FINS D'EVITER DES PROBLEMES DANS LA DECOMPOSITION
  673. if(iimpi.eq.1524) WRITE(IOIMP,7391) DMAX,DDDD,IENMU,
  674. & IENMU1,I,IEL
  675. 7391 FORMAT(' DMAX DDDD IENMU IENMU1 I IEL',2E12.5,4I3)
  676.  
  677. IF(DMAX.LE.xpetit.AND.IENMU.NE.IENMU1
  678. & .AND.IEL.EQ.1) GOTO 377
  679. IF(DMAX.LE.xzprec*dmaxge) DMAX = DMAXGE
  680. DMAX=DMAX*1.5D0
  681.  
  682. C XMATRI=IMATTT(IEL)
  683. C SEGACT,XMATRI
  684. C ... DMAXY = maximum des termes dans la première colonne (les 2 premiers exclus) ...
  685. DMAXY=sqrt(xpetit)*1d5
  686. if (norinc.eq.0) dmaxy=1.d0
  687. DO ICO=1,N1
  688. DMAXY = MAX ( DMAXY, ABS( RE(ICO,1,iel)))
  689. enddo
  690. DO iko=1,N3
  691. DMAXY = MAX ( DMAXY, ABS( RE(1,iko,iel)))
  692. enddo
  693.  
  694. DMAX = DMAX / DMAXY
  695. * Attention
  696. * on reprend le dmax du premier passage pour appliquer le meme facteur sur les lignes et les colonnes
  697. * car le max de la ligne n'est pas forcement celui de la colinne
  698. if (.not.ldiag) dmax=dnor(ival(1))
  699. ENDIF
  700. IF( IIMPI. EQ.1524 ) WRITE(IOIMP,7398) DMAX
  701. 7398 FORMAT(' facteur multiplicatif de norme ',e12.5)
  702.  
  703. C ... Copie de la matrice élémentaire fois DMAX*COER dans RA ...
  704. DO 21 ICO=1,N1
  705. DO 2110 IKO=1,N3
  706. RA(ICO,IKO)=RE(ICO,IKO,iel)*DMAX*COER
  707. 2110 CONTINUE
  708. 21 CONTINUE
  709. C ... La sous-matrice de dimension 2 est multipliée par DMAXY ...
  710. ** si on ne booste pas l'egalite des mults on a des problemes de precision sur ceux ci
  711. if (norinc.eq.0) dmaxy=dmaxy*2.d0
  712. RA(1,1)=RA(1,1)*DMAXY
  713. RA(2,1)=RA(2,1)*DMAXY
  714. RA(1,2)=RA(1,2)*DMAXY
  715. RA(2,2)=RA(2,2)*DMAXY
  716.  
  717. C ... Mise à DMAX des DNOR correspondant aux multiplicateurs ...
  718. IF(LDIAG) THEN
  719. DO 22 ICO=1,2
  720. DNOR(IVAL(ICO))=DMAX
  721. ** write (6,*) 'il dnor ',ival(ico),dnor(ival(ico))
  722. mdno1.DNOR(IVAL(ICO))=DMAX
  723. 22 CONTINUE
  724. ENDIF
  725. CMB ... Je suppose que l'on s'en servira pour multiplier les composantes
  726. C du second membre pour que les déplacements imposés soient corrects ...
  727.  
  728. C ... Boucle sur les variables primales ...
  729. ** write(6,*) 'ldmt2 trsup ',trsup
  730. DO 24 ICO=1,N3
  731. C ... INO = N° local du noeud ...
  732. INO=INUINV(NUM(NOELEP(ICO),IEL))
  733. C ... IO = N° de ligne du DDL N° ICO ...
  734. IO=IVAL(ICO)
  735. if (ico.eq.1) io1=io
  736. if (ico.eq.2) io2=io
  737. C ... LLIGN contient les lignes liées aux noeud INO ...
  738. LLIGN=ILIGN(INO)
  739. SEGACT,LLIGN*MOD
  740. C ... On rajoute le terme diagonal au DIAG ...
  741. IF(LDIAG) DIAG(IO)=DIAG(IO)+RA(ICO,ICO)
  742. C ... IMMMM contient les N°s des lignes des DDL ...
  743. DO 132 JLIJ=1,IMMMM(/1)
  744. JLIJ1=JLIJ
  745. C ... On est censé trouver le bon N° de ligne et aller au 133 ...
  746. IF( IMMMM(JLIJ).EQ.IO) GO TO 133
  747. 132 CONTINUE
  748.  
  749. IF(IIMPI.EQ.1524) WRITE(IOIMP,7354)
  750. 7354 FORMAT( ' PREMIERE ERREUR 5')
  751. CALL ERREUR(5)
  752. RETURN
  753.  
  754. 133 CONTINUE
  755.  
  756. C ... Deuxième niveau de boucle sur les variables primales ...
  757. DO 26 IRO=1,N1
  758. C ... IA = N° de colonne du DDL N° IRO ...
  759. IA=IVAL(IRO)
  760. C ... Cas symétrique !!! ...
  761. IF(IA.GT.IO) GO TO 26
  762.  
  763. C ... JLD et JLT sont égaux au début et à la fin (dans XXVA et LINC)
  764. C des données relatives au DDL N° IRO ...
  765. JLT=IPPO(JLIJ1+1)
  766. JLD=IPPO(JLIJ1)+1
  767.  
  768. C ... On cherche le bon N° de colonne dans LINC ...
  769. DO 134 JL=JLD,JLT
  770. JL1=JL
  771. IF(LINC(JL).EQ.IA) GO TO 135
  772. 134 CONTINUE
  773.  
  774. IF(IIMPI.EQ.1524) WRITE(IOIMP,7355)
  775. 7355 FORMAT( ' DEUXIEME ERREUR 5')
  776. CALL ERREUR(5)
  777. RETURN
  778.  
  779. 135 CONTINUE
  780. C ... On rajoute le coefficient dans RA au XXVA correspondant ...
  781. if(.not.trsup) then
  782. XXVA(JL1)=XXVA(JL1)+RA(ICO,IRO)
  783. else
  784. XXVA(JL1)=XXVA(JL1)+RA(IRO,ICO)
  785. endif
  786. 26 CONTINUE
  787. SEGDES,LLIGN
  788. 24 CONTINUE
  789. C on stocke dans ittr les couples de LX
  790. ittr(io1)=io2
  791. ittr(io2)=io1
  792. C SEGDES,XMATRI
  793. 14 CONTINUE
  794. C ... Après son assemblage, une matrice de blocage n'est plus stigmatisée comme telle ...
  795. INOMUL(I)=.TRUE.
  796. 377 CONTINUE
  797. SEGSUP,RA
  798. 11 CONTINUE
  799. C ... Fin de la boucle sur les rigidités ...
  800. GO TO 375
  801. C ... Fin de la boucle cachée ...
  802. 3750 CONTINUE
  803. ** do i = 1,inc
  804. ** if (ittr(i).eq.0.and.diag(i).ne.vmax(i))
  805. ** > write (6,*) 'diag vmax ',i,diag(i),vmax(i)
  806. ** enddo
  807.  
  808.  
  809.  
  810. C ... C'est l'ménaage finaale ...
  811. DO 18 IK=1,NNR
  812. INCTRA=INCTRR(IK)
  813. IF(MENAGE) THEN
  814. SEGSUP,INCTRA
  815. ELSE
  816. SEGDES,INCTRA
  817. ENDIF
  818. 18 CONTINUE
  819. C PV ON DESACTIVE TOUT
  820. NNR=IRIGEL(/2)
  821. DO 2 IRI=1,NNR
  822. DESCR=IRIGEL(3,IRI)
  823. SEGDES DESCR
  824. IPT1=IRIGEL(1,IRI)
  825. SEGDES IPT1
  826. xMATRI=IRIGEL(4,IRI)
  827. SEGDES xMATRI
  828. 2 CONTINUE
  829. INTERR(1)=NJTOT
  830. c-old IF(IPASG.EQ.0) CALL ERREUR(-278)
  831. IF(MENAGE.AND.IPASG.EQ.0) CALL ERREUR(-278)
  832. IPASG=1
  833.  
  834. IF(IIMPI.EQ.1457) WRITE(IOIMP,4821) LLVNUL,NJTOT
  835. 4821 FORMAT(' NB DE VALEURS NON NULLES DANS LA MATRICE ',I9,/
  836. & ' NB DE VALEURS DANS LA MATRICE ',I9)
  837. SEGSUP,IITOP,ITOPO
  838. IF(MENAGE) THEN
  839. C SEGSUP,IITOP
  840. ccc SEGSUP,IMINI
  841. SEGSUP,INUINV
  842. C SEGSUP,ITOPO
  843. SEGSUP,INCTRR,INCTRD
  844. ELSE
  845. C SEGDES,IITOP
  846. ccc SEGDES,IMINI
  847. SEGDES,INUINV
  848. C SEGDES,ITOPO
  849. SEGDES,INCTRR
  850. ENDIF
  851. C write(6,*) ' diag'
  852. C write(6,*) ( diag(iou),iou=1,diag(/1))
  853. SEGDES,MDIAG
  854. SEGDES,MIMIK
  855. IF (MDNOR.NE.0) SEGDES,MDNOR
  856. SEGDES,MILIGN
  857. SEGDES,MMATRI
  858. MMMTRI=MMATRI
  859. SEGDES,MINCPO
  860. SEGDES,MIPO1
  861. SEGDES,IPOS
  862. SEGSUP,IVAL,ITRA,JNOMUL
  863. IF(NORINC.NE.0 .AND. MENAGE) THEN
  864. CALL ASSEM0(MRIGID,2,inwuit)
  865. SEGDES,MRIGID
  866. ELSE
  867. SEGDES,MRIGID
  868. ENDIF
  869. RETURN
  870. END
  871.  
  872.  
  873.  
  874.  
  875.  
  876.  
  877.  
  878.  
  879.  
  880.  
  881.  
  882.  
  883.  
  884.  

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