Télécharger ldmt2.eso

Retour à la liste

Numérotation des lignes :

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

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