Télécharger assem5.eso

Retour à la liste

Numérotation des lignes :

  1. C ASSEM5 SOURCE PV 16/11/17 21:58:15 9180
  2. c assem5 source chat 94/08/12 21:15:05 1247
  3. SUBROUTINE ASSEM5(ITRAV1,ITOPO1,INUIN1,MMMTRI,INCTR1
  4. &,IITOP1,NBNNMA,SNTT)
  5. c
  6. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  7. c
  8. c **** subroutine pour faire l'assemblage de matrices symetriques
  9. c en vue d'une decomposition de crout modifiee.
  10. c
  11. c en entree:
  12. c **** itrav1 : pointeur objet mrigidite
  13. c **** itopo1 : pointeur segment de travail itopo ( voir assem1)
  14. c **** iitop1 : pointeur segment de travail iitop ( voir assem1)
  15. c **** inuin1 : pointeur segment de travail inuinv( voir assem1)
  16. c **** mmmtri : pointeur objet matrice triangularisee (non modifie)
  17. c (voir smmatri)
  18. c **** nbnnma : nombre d'inconnues maitresses
  19. c **** sntt : segment des noeuds non totalement maitres
  20. c
  21. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  22. c
  23. IMPLICIT REAL*8 (A-H,O-Z)
  24. IMPLICIT INTEGER(I-N)
  25. -INC CCREEL
  26. -INC CCOPTIO
  27. -INC SMELEME
  28. -INC SMRIGID
  29. -INC SMMATRI
  30. c
  31. SEGMENT,INUINV(NNGLOB)
  32. SEGMENT,ITOPO(IENNO)
  33. SEGMENT,IITOP(NNOE+1)
  34. SEGMENT,INCTRR(NIRI)
  35. SEGMENT,INCTRA(NLIGRE)
  36. SEGMENT,IPV(NNOE)
  37. SEGMENT,VMAX(INC)
  38. SEGMENT SNTT
  39. INTEGER NTTMAI(NN)
  40. ENDSEGMENT
  41. c
  42. c **** ces tableaux servent au reperage de la matrice pour l'assemblag
  43. c **** il seront tous supprimes en fin d'assemblage.
  44. c
  45. SEGMENT,IVAL(NNN)
  46. SEGMENT,ITRA(NNN,2)
  47. SEGMENT TRATRA
  48. REAL*8 XTRA(INCRED,INCDIF)
  49. INTEGER LTRA(INC,INCDIF)
  50. INTEGER NTRA(INCRED,INCDIF)
  51. INTEGER MTRA(INCDIF)
  52. ENDSEGMENT
  53. c
  54. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  55. c
  56. c **** ival(i)=j : la i eme ligne d'une petite matrice s'assemble
  57. c dans la j eme de la grande.
  58. c **** itrav(i,1)=j : la ieme inconnue du noeud en cours d'assemblage
  59. c et qui se trouve dans la petite matrice se trouve
  60. c en j eme position de la petite matrice.
  61. c **** itrav(i,2) : la ieme inconnue du noeud en cours d'assemblage
  62. c present dans la petite matrice est en jeme
  63. c position dans la grande
  64. c
  65. cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  66. c
  67. SEGMENT,RA(N1,N1)*D
  68. SEGMENT JNOMUL
  69. LOGICAL INOMUL(NNR)
  70. ENDSEGMENT
  71. REAL*8 DMAX,COER,DDDD,DMAXY
  72. LOGICAL NOMUL
  73. *
  74. * en cas de normalisation des variables.
  75. *
  76. INWUIT=0
  77. IF(NORINC.NE.0) CALL ASSEM0(ITRAV1,1,INWUIT)
  78. *
  79. * pv on active une fois pour toutes les meleme descr... de la rigidite
  80. * on en profite pour creer inomul
  81. c
  82. c **** recherche de la dimension max de ival,et segini de ival et itra
  83. c
  84. c
  85. INCTRR=INCTR1
  86. SEGACT,INCTRR
  87. c
  88. MRIGID=ITRAV1
  89.  
  90. SEGACT,MRIGID
  91. NNR=IRIGEL(/2)
  92. NNN=0
  93. SEGINI JNOMUL
  94. DO 1 IRI=1,NNR
  95. c
  96. INCTRA=INCTRR(IRI)
  97. SEGACT INCTRA
  98. c
  99. DESCR=IRIGEL(3,IRI)
  100. SEGACT, DESCR
  101. c
  102. IPT1=IRIGEL(1,IRI)
  103. SEGACT IPT1
  104. c
  105. XMATRI=IRIGEL(4,IRI)
  106. SEGACT XMATRI
  107. c
  108. c ... na = nb de variables primales ...
  109. NA=LISINC(/2)
  110. c ... nnn = max des nb de variables primales ...
  111. NNN=MAX(NA,NNN)
  112. c
  113. c ... inomul (non multiplicateurs ?) dit s'il ne s'agit pas des cl ...
  114. INOMUL(IRI)=.TRUE.
  115. IF(IPT1.ITYPEL.EQ.22) INOMUL(IRI)=.FALSE.
  116. c
  117. 1 CONTINUE
  118. c
  119. c ... ival a la taille nnn ...
  120. SEGINI,IVAL
  121. c ... itra a la taille (nnn,2) ...
  122. SEGINI,ITRA
  123. c
  124. c **** activation des segments de travails et de mmatri
  125. c
  126. c ... recherche du nombre d'inconnues impliquees
  127. c
  128. MMATRI=MMMTRI
  129. SEGACT,MMATRI*MOD
  130. c
  131. MINCPO=IINCPO
  132. SEGACT,MINCPO
  133. c
  134. c
  135. MAXIPO=0
  136. DO 150 I=1,INCPO(/2)
  137. DO 150 J=1,INCPO(/1)
  138. MAXIPO=MAX(INCPO(J,I),MAXIPO)
  139. 150 CONTINUE
  140. c
  141. N1=MAXIPO
  142. c
  143. ITOPO=ITOPO1
  144. SEGACT,ITOPO
  145. c
  146. IITOP=IITOP1
  147. SEGACT,IITOP
  148. c
  149. INUINV=INUIN1
  150. SEGACT,INUINV
  151. c
  152. c ... nnoe = nombre de noeuds impliqués ...
  153. c
  154. SEGACT SNTT
  155. NNOE=INCPO(/2)+NTTMAI(/1)
  156. c
  157. c ... inc = dimension de la matrice ...
  158. INC=MAXIPO
  159. c
  160. SEGINI,MDIAG
  161. IDIAG=MDIAG
  162. c
  163. SEGINI,MILIGN
  164. IILIGN=MILIGN
  165. c
  166. c ... incdif = nb de variables primales différentes ...
  167. INCDIF=INCPO(/1)
  168. c
  169. MIMIK=IIMIK
  170. SEGACT,MIMIK
  171. c
  172. NNOE=INCPO(/2)
  173. c
  174. c ... ipv est de taille nnoe ...
  175. SEGINI IPV
  176. c
  177. INCRED=0
  178. DO 80 INO =1,NNOE
  179. ICOMPT=0
  180. c ... maxele = nb d'éléments qui contiennent le noeud no ino ...
  181. MAXELE = (IITOP(INO+1)-IITOP(INO))/2
  182. c ... boucle sur ces éléments ...
  183. DO 81 IELE=1,MAXELE
  184. c
  185. c ... iiu = pointeur dans iitop ...
  186. IIU=IITOP(INO) + IELE + IELE -2
  187. c
  188. c ... iel = numero de l'élément dans son maillage ...
  189. IEL=ITOPO(IIU)
  190. c
  191. c ... iri = numero du maillage (dans irigel) de cet élément ...
  192. IRI=ITOPO(IIU+1)
  193. c
  194. c ... on va chercher ce maillage dans irigel ...
  195. MELEME=IRIGEL(1,IRI)
  196. c
  197. c ... boucle sur les noeuds de cet élément ...
  198. DO 83 I=1,NUM(/1)
  199. c
  200. c ... ip = numéro local du i-ème noeud de l'élément ...
  201. IP=INUINV(NUM(I,IEL))
  202. c
  203. c ... la première condition = les "connections" sont symétriques ...
  204. **pv IF (IP.GT.INO) GOTO 83
  205. c ... le noeud ip a déjà vu ce noeud ...
  206. IF (IPV(IP).EQ.INO) GOTO 83
  207. c
  208. c ... ces deux opérations se font si (ip <= ino et ipv(ip) != ino) ...
  209. c ... ipv = n° max des noeuds "connectés" au noeud ip ...
  210. c
  211. IPV(IP)=INO
  212. c ... icompt = nombre de noeuds connectés (de numéro < ou =) au noeud ino ...
  213. ICOMPT=ICOMPT+1
  214. 83 CONTINUE
  215. 81 CONTINUE
  216. c ... incred = maximum des icompt sur tous les noeuds impliqués ...
  217. c
  218. INCRED=MAX(INCRED,ICOMPT)
  219. 80 CONTINUE
  220. c
  221. SEGSUP IPV
  222. c
  223. INCRED=INCRED*INCDIF
  224. c ... la taille de tratra dépend de inc, incred et incdif ...
  225. SEGINI TRATRA
  226. segini vmax
  227. c
  228. c ... coefficients de normalisation ...
  229. c -----------------------------
  230. c
  231. LLVNUL=0
  232. IJMAX=0
  233. NJTOT=0
  234. c
  235. c ... la taille de mdnor est inc (nombre d'équations) ...
  236. SEGINI MDNOR
  237. IDNORM=MDNOR
  238. DO 20 I=1,INC
  239. DNOR(I)=1.D0
  240. 20 CONTINUE
  241. c
  242. c ... si une normalisation est imposée, les valeurs dans mdnor seront modifiées
  243. c par assem3 ...
  244. IF(NORINC.NE.0) CALL ASSEM3(MDNOR,MIMIK,MINCPO)
  245. c
  246. c **** boucle *100* sur les numeros de noeuds que l'on assemble
  247. c -------------------------------------------------------------
  248. c
  249. c ... ipre = n° de la première ligne concernée par le noeud ino ...
  250. IPRE=1
  251. c
  252. DO 999 IJK=1,2
  253. c
  254. DO 100 INO=1,NNOE
  255. DO 101 IIT=1,INCDIF
  256. MTRA(IIT)=0
  257. 101 CONTINUE
  258. c ... ider = n° de la dernière ligne concernée par le noeud ino ...
  259. IDER=0
  260. c
  261. IF(IJK.EQ.2) IPRE=MAXIPO+1
  262. DO 998 II=1,INCDIF
  263. IF (IJK.EQ.1.AND.INCPO(II,INO).LE.NBNNMA)
  264. # IDER=MAX(INCPO(II,INO),IDER)
  265. IF (IJK.EQ.2.AND.INCPO(II,INO).GT.NBNNMA) THEN
  266. IDER=MAX(INCPO(II,INO),IDER)
  267. IPRE=MIN(INCPO(II,INO),IPRE)
  268. ENDIF
  269. 998 CONTINUE
  270. IF (IDER.EQ.0.OR.IPRE.EQ.(MAXIPO+1)) GOTO 100
  271. LLVVA=0
  272. c
  273. c **** boucle *99* sur les elements touchant le noeud ino
  274. c pour les elements multiplicateur on ne fait pas
  275. c l'assemblage
  276. c
  277. c ... maxele = nb d'éléments qui contiennent le noeud no ino ...
  278. MAXELE= (IITOP(INO+1) -IITOP(INO))/2
  279. c
  280.  
  281. DO 99 IELE=1,MAXELE
  282. IIU=IITOP(INO) + IELE + IELE - 2
  283. c ... iel = numero de l'élément dans son maillage ...
  284. IEL=ITOPO(IIU)
  285. c ... iri = numero du maillage (dans irigel) de cet élément ...
  286. IRI=ITOPO(IIU+1)
  287. c ... meleme = pointeur vers ce maillage ...
  288. MELEME=IRIGEL(1,IRI)
  289. c ... descr = pointeur vers le segment descr concerné ...
  290. DESCR=IRIGEL(3,IRI)
  291. c ... inctra contient les indices (dans imik et ihar) des descriptions
  292. c des ddl correspondants ...
  293. INCTRA=INCTRR(IRI)
  294. c ... xmatri = tableau de pointeurs vers les matrices élémentaires ...
  295. XMATRI=IRIGEL(4,IRI)
  296.  
  297. c ... coer = coefficient multiplicateur ...
  298. COER=COERIG(IRI)
  299. c
  300. c **** nomul =.false. il existe un multiplicateur
  301. c **** initialisation de ival. ival(i)=j veut dire que
  302. c **** la i eme ligne de la petite matrice s'assemble dans
  303. c **** la j eme de la grande matrice.
  304. c
  305. c ... nin = nombre de variables primales (dans descr) ...
  306. NIN=LISINC(/2)
  307. c ... non multiplicateur ? ...
  308. NOMUL=INOMUL(IRI)
  309. NA=0
  310. c ... boucle sur les variables primales de l'élément ...
  311. c
  312. DO 98 ICO=1,NIN
  313. c ... ija = numéro local du noeud-support du ddl ...
  314. IJA=INUINV(NUM(NOELEP(ICO),IEL))
  315. c ... ijb = indice du ddl ...
  316. IJB=INCTRA(ICO)
  317. c ... on met dans ival(ico) le n° de l'équation correspondant au noeud ija
  318. c et ddl no ijb ...
  319. IVAL(ICO)=INCPO(IJB,IJA)
  320. c ... si on a trouvé le bon noeud ...
  321. IF(IJA.EQ.INO) THEN
  322. c ... on incrémente na (le nombre de ddl connectés au noeud ino ?) ...
  323. NA=NA+1
  324. ITRA(NA,1)=ICO
  325. ITRA(NA,2)=IVAL(ICO)
  326. ENDIF
  327. 98 CONTINUE
  328. * XMATRI=IMATTT(IEL)
  329. * SEGACT,XMATRI
  330. c
  331. c **** boucle *85* sur les inconnues de la petite matrice
  332. c
  333. DO 95 IK=1,NIN
  334. c ... io = n° de l'équation correspondant au ddl no ik de l'élément (n° de la colonne) ...
  335. IO=IVAL(IK)
  336. c ... s'il dépasse le dernier intéressant on ne fait rien ...
  337. IF(IO.GT.IDER) GO TO 95
  338. c ... boucle sur les ddl's du noeud ino ...
  339. DO 90 INCC=1,NA
  340. c ... inco = le n° d'équation (ou de ligne) du ddl no incc du noeud ino ...
  341. INCO=ITRA(INCC,2)
  342. c
  343. IF (IJK.EQ.1) THEN
  344. IF (INCO.GT.IDER) GOTO 90
  345. ENDIF
  346. c
  347. c ... si io > inco on ne fait plus rien dans cette boucle ...
  348. c ... car c'est le cas symétrique !!! ...
  349. IF (IO.GT.INCO) GO TO 90
  350. c
  351. c ... iloc = n° de ligne relatif / ipre du noeud ino ...
  352. c
  353. ILOC=INCO-IPRE+1
  354. IF (IJK.EQ.2) THEN
  355. IF (ILOC.LE.0) GOTO 90
  356. ENDIF
  357. c
  358. c ... jj = n° (dans l'élément) de la colonne (et de la ligne) correspondant au
  359. c ddl n° incc du noeud ino ...
  360. JJ=ITRA(INCC,1)
  361. c ... ci-dessous on stockera les quelques lignes de la rigidité concernant le
  362. c noeud ino dans le segment tratra. les coefficients vont dans xtra, le deuxième
  363. c indice (iloc) est le numéro relatif de la ligne, le premier (imtt) n'a rien à voir avec le
  364. c numéro de la colonne. celui-ci est stocké dans ntra(imtt,iloc). l'information inverse
  365. c (imtt en fonction de io) se trouve dans ltra(io,iloc). on met dans mtra le nombre
  366. c de termes différents assemblés par ligne ...
  367. ILTT= LTRA(IO,ILOC)
  368. IF(ILTT.EQ.0) THEN
  369. c ... llvva est initialisé à 0 au début de la boucle 100 ...
  370. LLVVA=LLVVA+1
  371. c ... les mtra sont mis à 0 au début de la boucle 100 ...
  372. IMMTT=MTRA(ILOC)+1
  373. MTRA(ILOC)=IMMTT
  374. XTRA(IMMTT,ILOC)=0.D0
  375. NTRA(IMMTT,ILOC)=IO
  376. LTRA(IO,ILOC)=IMMTT
  377. ILTT=IMMTT
  378. ENDIF
  379. c
  380. c ... si ce n'est pas une condition aux limites, on assemble ...
  381. IF(NOMUL) THEN
  382. cfaux !!! xtra(iltt,iloc)=xtra(iltt,iloc)+re(ik,jj)*coer
  383. c ... cette ligne a été remplacée par ce qui suit car pour les termes
  384. c de couplage entre les ddl du même noeud, les termes au-dessus de
  385. c la diagonale étaient pris. Ça marchait car les matrices étaient
  386. c symétriques, comme maintenant on travaillera aussi sur des asymétriques ...
  387. IF(JJ.LT.IK) THEN
  388. IIILIG=IK
  389. IIICOL=JJ
  390. ELSE
  391. IIILIG=JJ
  392. IIICOL=IK
  393. ENDIF
  394. c ... attention ! dans xtra 1er indice est lié à la colonne, deuxième est le numéro
  395. c de la ligne, tandis que dans re c'est l'inverse ...
  396. XTRA(ILTT,ILOC)=XTRA(ILTT,ILOC)+RE(IIILIG,IIICOL,IEL)
  397. & *COER
  398. ENDIF
  399. 90 CONTINUE
  400. 95 CONTINUE
  401. * SEGDES,XMATRI
  402. 99 CONTINUE
  403. c
  404. c *** compactage des lignes, en meme temps calcul de ijmax qui sera
  405. c *** la dimension max d'un segment lign.
  406. c *** le segment associe a une ligne (segment llign) est de la forme :
  407. c *** immmm(na) permet de savoir si un mouvement d'ensemble sur la
  408. c *** ligne existe. ippo(na+1) donne la position dans xxva la 1ere
  409. c *** valeur de la ligne .xxva valeur de la matrice.
  410. c *** linc(i) donne le numero de la colonne du ieme elem de xxva
  411. c
  412. c ... na = nb de lignes concernant le noeud ino ...
  413. NA = IDER-IPRE+1
  414. c ... llvnul = somme cumulée des llvva ...
  415. LLVNUL=LLVNUL+LLVVA
  416. c
  417. SEGINI,LLIGN
  418. c
  419. NBINO=INO
  420. IF (IJK.EQ.2) THEN
  421. DO 116 JJ=1,NTTMAI(/1)
  422. IF (NTTMAI(JJ).EQ.INO) GOTO 117
  423. 116 CONTINUE
  424. GOTO 118
  425. 117 CONTINUE
  426. NBINO=NNOE+JJ
  427. ENDIF
  428. 118 CONTINUE
  429. ILIGN(NBINO)=LLIGN
  430. NBA=0
  431. c
  432. c ... boucle sur les ddl du noeud ino ...
  433. DO 120 JPA=1,NA
  434. c ... iiin = n° de ligne du ddl n° jpa ...
  435. IIIN=IPRE+JPA -1
  436. c
  437. c ... que l'on met au bon endroit dans immmm (faisant partie de llign) ...
  438. c
  439. IMMMM(JPA)=IIIN
  440. c
  441. c
  442. c ... ippo(i)+1 = adresse du début (dans xxva et linc) des données relatives au ddl n° i ...
  443. IPPO(JPA)=NBA
  444. c
  445. c ... boucle sur les termes dans la ligne ...
  446. DO 121 IPAK = 1,MTRA(JPA)
  447. c ... iunpak = n° de la colonne du terme n° ipak ...
  448. IUNPAK=NTRA(IPAK,JPA)
  449. c ... pour les termes mis dans llign on efface l'information sur la position dans xtra ...
  450. LTRA(IUNPAK,JPA)=0
  451. c ... compteur ++ ...
  452. NBA=NBA+1
  453. c ... le n° de la colonne va dans linc ...
  454. LINC(NBA)=IUNPAK
  455. c ... transfert du xtra (segment tratra) vers xxva (segment llign) ...
  456. XXVA(NBA)=XTRA(IPAK,JPA)
  457. vmax(iiin)=max(abs(xxva(nba)),vmax(iiin))
  458. c ... les termes diagonaux vont dans diag (segment mdiag) ...
  459. IF(IIIN.EQ.IUNPAK) DIAG(IIIN)=XXVA(NBA)
  460. 121 CONTINUE
  461. 120 CONTINUE
  462. c ... le pointeur vers la fin de la dernière ligne ...
  463. IPPO(NA+1)= NBA
  464. NJMAX= 0
  465. * recherche du mini globale sur toutes les inconnues
  466. LPA=IPRE
  467. c ... boucle sur les lignes relatives au noeud ino ...
  468. DO 126 JPA=IPRE,IDER
  469. c ... on met le n° du noeud dans ipno (segment milign) ...
  470. IPNO(JPA)=NBINO
  471. c ... ipde et ipdf : début et fin des données relatives au ddl n° jpa dans xxva et linc ...
  472. IPDE=IPPO(JPA-IPRE+1)+1
  473. IPDF=IPPO(JPA-IPRE+2)
  474. c ... lpa = n° mini de la colonne avec des termes non nuls dans la ligne n° jpa ...
  475. DO 155 JHT=IPDE,IPDF
  476. LPA=MIN(LPA,LINC(JHT))
  477. 155 CONTINUE
  478. c ... on le met dans ldeb (segment llign) ...
  479. LDEB(JPA-IPRE+1)=LPA
  480. c ... nna = longueur de la "skyline" ...
  481. NNA= JPA- LPA +1
  482. c ... njmax = profil cumulé sur un noeud ...
  483. NJMAX=NJMAX+NNA
  484. 126 CONTINUE
  485. DO 127 JPA=IPRE,IDER
  486. LDEB(JPA-IPRE+1)=LPA
  487. NNA= JPA- LPA +1
  488. NJMAX=NJMAX+NNA
  489. 127 continue
  490. c
  491. c ... njtot = profil total ...
  492. NJTOT=NJTOT+NJMAX
  493. c ... ijmax = (profil / noeud) maxi ...
  494. IF(IJMAX.LT.NJMAX) IJMAX=NJMAX
  495. SEGDES,LLIGN
  496. IPRE=IDER+1
  497. 100 CONTINUE
  498. c
  499. 999 CONTINUE
  500. c
  501. c ... fin de la boucle sur les noeuds ...
  502. SEGSUP TRATRA
  503. C
  504. C **** ON REPREND TOUTE LES MATRICES CONTENANT LES MULTIPLICATEURS
  505. C **** POUR MULTIPLIER TOUS LEURS TERMES PAR UNE NORME ATTACHEE
  506. C **** A CHAQUE MULTIPLICATEUR.PUIS ON LES ASSEMBLE.
  507. C
  508. * d'abord etablir une norme generale pour le cas ou on n'arrive pas
  509. * a calculer la norme particuliere
  510. DMAXGE=xpetit
  511. DO 378 I=1,INC
  512. ** write (6,*) ' assem5 diag vmav ',diag(i),vmax(i)
  513. DMAXGE=MAX(DMAXGE,vmax(i))
  514. 378 CONTINUE
  515. if (dmaxge*xzprec.lt.xpetit*10) dmaxge=1.d0
  516. if (iimpi.eq.1524)
  517. > write (6,*) ' nb inconnues facteur multiplicatif general ',
  518. > INC,DMAXGE
  519. IENMU=0
  520. c ... attention ! ici commence une boucle cachée (exécutée avec un goto 375) ...
  521. 375 IENMU1 = IENMU
  522. IENMU=0
  523. c ... boucle sur les rigidités qui calcule le nombre de matrices de blocages ...
  524. DO 376 I=1,NNR
  525. IF(.NOT.INOMUL(I)) IENMU=IENMU+1
  526. 376 CONTINUE
  527. c
  528. c ... s'il n'y en a pas on saute cette partie du code ...
  529. IF( IENMU.EQ.0) GO TO 3750
  530. c
  531. c ... mimik contient les noms des variables primales ...
  532. MIMIK=IIMIK
  533. SEGACT,MIMIK
  534. c ... boucle sur les rigidités ...
  535. DO 11 I=1,NNR
  536. c ... si celle si n'est pas une matrice de blocage on passe à la suivante ...
  537. IF(INOMUL(I)) GO TO 11
  538. c
  539. DESCR=IRIGEL(3,I)
  540. c ... n3 = nb de variables primales ...
  541. N3=LISINC(/2)
  542. COER=COERIG(I)
  543. MELEME=IRIGEL(1,I)
  544. INCTRA=INCTRR(I)
  545. XMATRI=IRIGEL(4,I)
  546. c ... n2 = nombre d'éléments dans le support géométrique ...
  547. N2=NUM(/2)
  548. c ... À quoi correspond ce cas ? (pas de matrices élémentaires) ...
  549. IF (RE(/3).EQ.0) THEN
  550. INOMUL(I)=.TRUE.
  551. GOTO 11
  552. ENDIF
  553. * XMATRI=IMATTT(1)
  554. * SEGACT,XMATRI
  555. c ... n1 = nombre de variables duales ...
  556. c ... pourquoi va-t-on chercher n3 dans descr et n1 dans re ? ...
  557. N1=RE(/1)
  558. SEGINI,RA
  559. c ... boucle sur les éléments ...
  560. DO 14 IEL=1,N2
  561. c ... boucle sur les inconnues ...
  562. DO 15 ICO=1,N3
  563. c ... ija = n° local du noeud-support de l'inconnue n° ico ...
  564. IJA=INUINV(NUM(NOELEP(ICO),IEL))
  565. c ... ijb = n° de l'inconnue ...
  566. IJB=INCTRA(ICO)
  567. c ... ival contient la correspondance entre le n° local du ddl
  568. c et le n° d'équation (de ligne) correspondant ...
  569. IVAL(ICO)=INCPO(IJB,IJA)
  570. 15 CONTINUE
  571.  
  572. DDDD= xpetit
  573. c ... dmax = max des valeurs absolues des termes diagonaux correspondants
  574. c au ddl de l'élément n° iel ...
  575. DMAX=DDDD
  576. if (n3.lt.3) then
  577. do ico=1,n3
  578. DMAX=MAX(DMAX,vmax(IVAL(ICO)))
  579. enddo
  580. endif
  581. DO 19 ICO=3,N3
  582. DMAX=MAX(DMAX,vmax(IVAL(ICO)))
  583. 19 CONTINUE
  584. ** write (6,*) ' assem2 dmax dmaxge ',dmax,dmaxge
  585. C AUX FINS D'EVITER DES PROBLEMES DANS LA DECOMPOSITION
  586. IF( IIMPI. EQ.1524 ) WRITE(IOIMP,7391)DMAX,DDDD,IENMU,IENMU1
  587. 1,I,IEL
  588. 7391 FORMAT(' DMAX DDDD IENMU IENMU1 I IEL',2E12.5,4I3)
  589. ** write (6,*) ' assem2 dmax dmxge',dmax,dmaxge
  590. IF(DMAX.LE.xzprec*dmaxge.AND.IENMU.NE.IENMU1.AND.IEL.EQ.1)GOTO 377
  591. IF(DMAX.LE.xzprec*dmaxge) DMAX = DMAXGE
  592. * facteur de normalisation cf PV pour ne pas avoir de pivot nul
  593. DMAX=DMAX*1.5
  594. * on penalise la matrice en cas de resolution iterative
  595. **pv if (nucrou.eq.1) DMAX=DMAX*1D5
  596. ** write (6,*) ' assem2 i iel dmax ',i,iel,dmax
  597. * XMATRI=IMATTT(IEL)
  598. * SEGACT,XMATRI
  599. c ... dmaxy = maximum des termes dans la première colonne (les 2 premiers exclus) ...
  600. DMAXY=SQRT(XPETIT)*1D5
  601. DO 821 ICO=3,N1
  602. DMAXY = MAX ( DMAXY, ABS( RE(ICO,1,iel)))
  603. 821 CONTINUE
  604. c
  605. IF( IIMPI. EQ.9022 ) WRITE(IOIMP,7398) DMAX
  606. 7398 FORMAT('FACTEUR MULTIPLICATIF DE NORME=',E12.5)
  607. c
  608. c ... copie de la matrice élémentaire fois dmax*coer dans ra ...
  609. DO 21 ICO=1,N1
  610. DO 2110 IKO=1,N1
  611. RA(ICO,IKO)=RE(ICO,IKO,iel)*DMAX*COER
  612. 2110 CONTINUE
  613. 21 CONTINUE
  614. c ... la sous-matrice de dimension 2 est multipliée par dmaxy ...
  615. * commenté car pose un problème dans excite et ne sert sans doute à rien
  616. * write (6,*) ' dmaxy dans assem5 ', dmaxy
  617. xcoef= 1D0
  618. RA(1,1)=RA(1,1)*xcoef
  619. RA(2,1)=RA(2,1)*xcoef
  620. RA(1,2)=RA(1,2)*xcoef
  621. RA(2,2)=RA(2,2)*xcoef
  622. c ... mise à dmax des dnor correspondant aux multiplicateurs ...
  623. DO 22 ICO=1,2
  624. c-old dnor(ival(ico))=dmax
  625. DNOR(IVAL(ICO))=DNOR(IVAL(ICO))*DMAX
  626. 22 CONTINUE
  627. c
  628. c ... boucle sur les variables primales ...
  629. DO 24 ICO=1,N3
  630. c ... ino = n° local du noeud ...
  631. INO=INUINV(NUM(NOELEP(ICO),IEL))
  632. c ... io = n° de ligne du ddl n° ico ...
  633. IO=IVAL(ICO)
  634. if (ico.eq.1) io1=io
  635. if (ico.eq.2) io2=io
  636. c ... cas des noeuds non totalement maitres ...
  637. IF(IO.GT.NBNNMA) THEN
  638. DO 1161 JJ=1,NTTMAI(/1)
  639. IF (NTTMAI(JJ).EQ.INO) GOTO 1171
  640. 1161 CONTINUE
  641. GOTO 1181
  642. 1171 CONTINUE
  643. INO=NNOE+JJ
  644. 1181 CONTINUE
  645. ENDIF
  646. c ... llign contient les lignes liées aux noeud ino ...
  647. LLIGN=ILIGN(INO)
  648. SEGACT,LLIGN*MOD
  649. c ... on rajoute le terme diagonal au diag ...
  650. DIAG(IO)=DIAG(IO)+RA(ICO,ICO)
  651. c ... immmm contient les n°s des lignes des ddl ...
  652. DO 132 JLIJ=1,IMMMM(/1)
  653. JLIJ1=JLIJ
  654. c ... on est censé trouver le bon n° de ligne et aller au 133 ...
  655. IF( IMMMM(JLIJ).EQ.IO) GO TO 133
  656. 132 CONTINUE
  657. c
  658. IF(IIMPI.EQ.1524) WRITE(IOIMP,7354)
  659. 7354 FORMAT( ' PREMIERE ERREUR 5')
  660. CALL ERREUR(5)
  661. RETURN
  662. c
  663. 133 CONTINUE
  664. c
  665. c ... deuxième niveau de boucle sur les variables primales ...
  666. DO 26 IRO=1,N3
  667. c ... ia = n° de colonne du ddl n° iro ...
  668. IA=IVAL(IRO)
  669. c ... cas symétrique !!! ...
  670. IF(IA.GT.IO) GO TO 26
  671. c
  672. c ... jld et jlt sont égaux au début et à la fin (dans xxva et linc)
  673. c des données relatives au ddl n° iro ...
  674. JLT=IPPO(JLIJ1+1)
  675. JLD=IPPO(JLIJ1)+1
  676. c
  677. c ... on cherche le bon n° de colonne dans linc ...
  678. DO 134 JL=JLD,JLT
  679. JL1=JL
  680. IF(LINC(JL).EQ.IA) GO TO 135
  681. 134 CONTINUE
  682. c
  683. IF(IIMPI.EQ.1524) WRITE(IOIMP,7355)
  684. 7355 FORMAT( ' DEUXIEME ERREUR 5')
  685. CALL ERREUR(5)
  686. RETURN
  687. c
  688. 135 CONTINUE
  689. c ... on rajoute le coefficient dans ra au xxva correspondant ...
  690. XXVA(JL1)=XXVA(JL1)+RA(ICO,IRO)
  691. 26 CONTINUE
  692. SEGDES,LLIGN
  693. 24 CONTINUE
  694. * on stocke dans ittr les couples de LX
  695. ittr(io1)=io2
  696. ittr(io2)=io1
  697. * SEGDES,XMATRI
  698. 14 CONTINUE
  699. c ... après son assemblage, une matrice de blocage n'est plus stigmatisée comme telle ...
  700. INOMUL(I)=.TRUE.
  701. 377 CONTINUE
  702. SEGSUP,RA
  703. 11 CONTINUE
  704. c ... fin de la boucle sur les rigidités ...
  705. GO TO 375
  706. c ... fin de la boucle cachée ...
  707. 3750 CONTINUE
  708. c
  709. DO 18 IK=1,NNR
  710. INCTRA=INCTRR(IK)
  711. SEGSUP,INCTRA
  712. 18 CONTINUE
  713. * pv on desactive tout
  714. NNR=IRIGEL(/2)
  715. DO 2 IRI=1,NNR
  716. DESCR=IRIGEL(3,IRI)
  717. SEGDES DESCR
  718. IPT1=IRIGEL(1,IRI)
  719. SEGDES IPT1
  720. xMATRI=IRIGEL(4,IRI)
  721. SEGDES xMATRI
  722. 2 CONTINUE
  723. INTERR(1)=NJTOT
  724. c
  725. IF(IIMPI.GE.1) WRITE(IOIMP,4821) LLVNUL,NJTOT
  726. 4821 FORMAT('Nbre de valeurs non nulles à l assemblage ',I9,/
  727. & 'Nbre de valeurs à l assemblage ',I9)
  728. c
  729. c
  730. SEGDES,SNTT
  731. SEGSUP,INCTRR
  732. SEGDES,MDIAG
  733. SEGDES,MIMIK
  734. SEGDES,MDNOR
  735. SEGSUP,ITOPO
  736. SEGSUP,IITOP
  737. SEGDES,MILIGN
  738. SEGSUP,INUINV
  739. SEGDES,MMATRI
  740. MMMTRI=MMATRI
  741. SEGDES,MINCPO
  742. SEGSUP,IVAL,ITRA,JNOMUL,vmax
  743. IF(NORINC.NE.0) THEN
  744. CALL ASSEM0(MRIGID,2,INWUIT)
  745. SEGDES,MRIGID
  746. ELSE
  747. SEGDES,MRIGID
  748. ENDIF
  749. RETURN
  750. END
  751. c
  752. c
  753. c
  754. c
  755. c
  756.  
  757.  
  758.  
  759.  
  760.  
  761.  
  762.  
  763.  
  764.  
  765.  
  766.  
  767.  
  768.  
  769.  
  770.  
  771.  
  772.  
  773.  
  774.  
  775.  
  776.  
  777.  

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