Télécharger assem5.eso

Retour à la liste

Numérotation des lignes :

  1. C ASSEM5 SOURCE PV 18/11/26 21:15:04 1009
  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 *95* sur les inconnues de la petite matrice
  332. c
  333. DO 90 INCC=1,NA
  334. c ... inco = le n° d'équation (ou de ligne) du ddl no incc du noeud ino ...
  335. INCO=ITRA(INCC,2)
  336. IF (IJK.EQ.1) then
  337. IF (INCO.GT.IDER) GOTO 90
  338. ENDIF
  339. c ... iloc = n° de ligne relatif / ipre du noeud ino ...
  340. c
  341. ILOC=INCO-IPRE+1
  342. IF (IJK.EQ.2) THEN
  343. IF (ILOC.LE.0) GOTO 90
  344. ENDIF
  345. c
  346. c ... jj = n° (dans l'élément) de la colonne (et de la ligne) correspondant au
  347. c ddl n° incc du noeud ino ...
  348. JJ=ITRA(INCC,1)
  349. DO 95 IK=1,NIN
  350. c ... io = n° de l'équation correspondant au ddl no ik de l'élément (n° de la colonne) ...
  351. IO=IVAL(IK)
  352. c ... s'il dépasse le dernier intéressant on ne fait rien ...
  353. IF(IO.GT.IDER) GO TO 95
  354. c ... boucle sur les ddl's du noeud ino ...
  355. c
  356. c
  357. c ... si io > inco on ne fait plus rien dans cette boucle ...
  358. c ... car c'est le cas symétrique !!! ...
  359. IF (IO.GT.INCO) GO TO 95
  360. c
  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.  
  397. ** XTRA(ILTT,ILOC)=XTRA(ILTT,ILOC)+RE(IIILIG,IIICOL,IEL)
  398. *** on exploite la symetrie de re
  399. XTRA(ILTT,ILOC)=XTRA(ILTT,ILOC)+RE(ik,jj,IEL)
  400. & *COER
  401. ENDIF
  402. 95 CONTINUE
  403. 90 CONTINUE
  404. * SEGDES,XMATRI
  405. 99 CONTINUE
  406. c
  407. c *** compactage des lignes, en meme temps calcul de ijmax qui sera
  408. c *** la dimension max d'un segment lign.
  409. c *** le segment associe a une ligne (segment llign) est de la forme :
  410. c *** immmm(na) permet de savoir si un mouvement d'ensemble sur la
  411. c *** ligne existe. ippo(na+1) donne la position dans xxva la 1ere
  412. c *** valeur de la ligne .xxva valeur de la matrice.
  413. c *** linc(i) donne le numero de la colonne du ieme elem de xxva
  414. c
  415. c ... na = nb de lignes concernant le noeud ino ...
  416. NA = IDER-IPRE+1
  417. c ... llvnul = somme cumulée des llvva ...
  418. LLVNUL=LLVNUL+LLVVA
  419. c
  420. SEGINI,LLIGN
  421. c
  422. NBINO=INO
  423. IF (IJK.EQ.2) THEN
  424. DO 116 JJ=1,NTTMAI(/1)
  425. IF (NTTMAI(JJ).EQ.INO) GOTO 117
  426. 116 CONTINUE
  427. GOTO 118
  428. 117 CONTINUE
  429. NBINO=NNOE+JJ
  430. ENDIF
  431. 118 CONTINUE
  432. ILIGN(NBINO)=LLIGN
  433. NBA=0
  434. c
  435. c ... boucle sur les ddl du noeud ino ...
  436. DO 120 JPA=1,NA
  437. c ... iiin = n° de ligne du ddl n° jpa ...
  438. IIIN=IPRE+JPA -1
  439. c
  440. c ... que l'on met au bon endroit dans immmm (faisant partie de llign) ...
  441. c
  442. IMMMM(JPA)=IIIN
  443. c
  444. c
  445. c ... ippo(i)+1 = adresse du début (dans xxva et linc) des données relatives au ddl n° i ...
  446. IPPO(JPA)=NBA
  447. c
  448. c ... boucle sur les termes dans la ligne ...
  449. DO 121 IPAK = 1,MTRA(JPA)
  450. c ... iunpak = n° de la colonne du terme n° ipak ...
  451. IUNPAK=NTRA(IPAK,JPA)
  452. c ... pour les termes mis dans llign on efface l'information sur la position dans xtra ...
  453. LTRA(IUNPAK,JPA)=0
  454. c ... compteur ++ ...
  455. NBA=NBA+1
  456. c ... le n° de la colonne va dans linc ...
  457. LINC(NBA)=IUNPAK
  458. c ... transfert du xtra (segment tratra) vers xxva (segment llign) ...
  459. XXVA(NBA)=XTRA(IPAK,JPA)
  460. vmax(iiin)=max(abs(xxva(nba)),vmax(iiin))
  461. c ... les termes diagonaux vont dans diag (segment mdiag) ...
  462. IF(IIIN.EQ.IUNPAK) DIAG(IIIN)=XXVA(NBA)
  463. 121 CONTINUE
  464. 120 CONTINUE
  465. c ... le pointeur vers la fin de la dernière ligne ...
  466. IPPO(NA+1)= NBA
  467. NJMAX= 0
  468. * recherche du mini globale sur toutes les inconnues
  469. LPA=IPRE
  470. c ... boucle sur les lignes relatives au noeud ino ...
  471. DO 126 JPA=IPRE,IDER
  472. c ... on met le n° du noeud dans ipno (segment milign) ...
  473. IPNO(JPA)=NBINO
  474. c ... ipde et ipdf : début et fin des données relatives au ddl n° jpa dans xxva et linc ...
  475. IPDE=IPPO(JPA-IPRE+1)+1
  476. IPDF=IPPO(JPA-IPRE+2)
  477. c ... lpa = n° mini de la colonne avec des termes non nuls dans la ligne n° jpa ...
  478. DO 155 JHT=IPDE,IPDF
  479. LPA=MIN(LPA,LINC(JHT))
  480. 155 CONTINUE
  481. c ... on le met dans ldeb (segment llign) ...
  482. LDEB(JPA-IPRE+1)=LPA
  483. c ... nna = longueur de la "skyline" ...
  484. NNA= JPA- LPA +1
  485. c ... njmax = profil cumulé sur un noeud ...
  486. NJMAX=NJMAX+NNA
  487. 126 CONTINUE
  488. DO 127 JPA=IPRE,IDER
  489. LDEB(JPA-IPRE+1)=LPA
  490. NNA= JPA- LPA +1
  491. NJMAX=NJMAX+NNA
  492. 127 continue
  493. c
  494. c ... njtot = profil total ...
  495. NJTOT=NJTOT+NJMAX
  496. c ... ijmax = (profil / noeud) maxi ...
  497. IF(IJMAX.LT.NJMAX) IJMAX=NJMAX
  498. SEGDES,LLIGN
  499. IPRE=IDER+1
  500. 100 CONTINUE
  501. c
  502. 999 CONTINUE
  503. c
  504. c ... fin de la boucle sur les noeuds ...
  505. SEGSUP TRATRA
  506. C
  507. C **** ON REPREND TOUTE LES MATRICES CONTENANT LES MULTIPLICATEURS
  508. C **** POUR MULTIPLIER TOUS LEURS TERMES PAR UNE NORME ATTACHEE
  509. C **** A CHAQUE MULTIPLICATEUR.PUIS ON LES ASSEMBLE.
  510. C
  511. * d'abord etablir une norme generale pour le cas ou on n'arrive pas
  512. * a calculer la norme particuliere
  513. DMAXGE=xpetit
  514. DO 378 I=1,INC
  515. ** write (6,*) ' assem5 diag vmav ',diag(i),vmax(i)
  516. DMAXGE=MAX(DMAXGE,abs(vmax(i)))
  517. 378 CONTINUE
  518. if (dmaxge.lt.xpetit/xzprec*10) dmaxge=1.d0
  519. if (iimpi.eq.1524)
  520. > write (6,*) ' nb inconnues facteur multiplicatif general ',
  521. > INC,DMAXGE
  522. IENMU=0
  523. c ... attention ! ici commence une boucle cachée (exécutée avec un goto 375) ...
  524. 375 IENMU1 = IENMU
  525. IENMU=0
  526. c ... boucle sur les rigidités qui calcule le nombre de matrices de blocages ...
  527. DO 376 I=1,NNR
  528. IF(.NOT.INOMUL(I)) IENMU=IENMU+1
  529. 376 CONTINUE
  530. c
  531. c ... s'il n'y en a pas on saute cette partie du code ...
  532. IF( IENMU.EQ.0) GO TO 3750
  533. c
  534. c ... mimik contient les noms des variables primales ...
  535. MIMIK=IIMIK
  536. SEGACT,MIMIK
  537. c ... boucle sur les rigidités ...
  538. DO 11 I=1,NNR
  539. c ... si celle si n'est pas une matrice de blocage on passe à la suivante ...
  540. IF(INOMUL(I)) GO TO 11
  541. c
  542. DESCR=IRIGEL(3,I)
  543. c ... n3 = nb de variables primales ...
  544. N3=LISINC(/2)
  545. COER=COERIG(I)
  546. MELEME=IRIGEL(1,I)
  547. INCTRA=INCTRR(I)
  548. XMATRI=IRIGEL(4,I)
  549. c ... n2 = nombre d'éléments dans le support géométrique ...
  550. N2=NUM(/2)
  551. c ... À quoi correspond ce cas ? (pas de matrices élémentaires) ...
  552. IF (RE(/3).EQ.0) THEN
  553. INOMUL(I)=.TRUE.
  554. GOTO 11
  555. ENDIF
  556. * XMATRI=IMATTT(1)
  557. * SEGACT,XMATRI
  558. c ... n1 = nombre de variables duales ...
  559. c ... pourquoi va-t-on chercher n3 dans descr et n1 dans re ? ...
  560. N1=RE(/1)
  561. SEGINI,RA
  562. c ... boucle sur les éléments ...
  563. DO 14 IEL=1,N2
  564. c ... boucle sur les inconnues ...
  565. DO 15 ICO=1,N3
  566. c ... ija = n° local du noeud-support de l'inconnue n° ico ...
  567. IJA=INUINV(NUM(NOELEP(ICO),IEL))
  568. c ... ijb = n° de l'inconnue ...
  569. IJB=INCTRA(ICO)
  570. c ... ival contient la correspondance entre le n° local du ddl
  571. c et le n° d'équation (de ligne) correspondant ...
  572. IVAL(ICO)=INCPO(IJB,IJA)
  573. 15 CONTINUE
  574.  
  575. DDDD= xpetit
  576. c ... dmax = max des valeurs absolues des termes diagonaux correspondants
  577. c au ddl de l'élément n° iel ...
  578. DMAX=DDDD
  579. if (n3.lt.3) then
  580. do ico=1,n3
  581. DMAX=MAX(DMAX,vmax(IVAL(ICO)))
  582. enddo
  583. endif
  584. DO 19 ICO=3,N3
  585. DMAX=MAX(DMAX,vmax(IVAL(ICO)))
  586. 19 CONTINUE
  587. ** write (6,*) ' assem2 dmax dmaxge ',dmax,dmaxge
  588. C AUX FINS D'EVITER DES PROBLEMES DANS LA DECOMPOSITION
  589. IF( IIMPI. EQ.1524 ) WRITE(IOIMP,7391)DMAX,DDDD,IENMU,IENMU1
  590. 1,I,IEL
  591. 7391 FORMAT(' DMAX DDDD IENMU IENMU1 I IEL',2E12.5,4I3)
  592. ** write (6,*) ' assem2 dmax dmxge',dmax,dmaxge
  593. IF(DMAX.LE.xzprec*dmaxge.AND.IENMU.NE.IENMU1.AND.IEL.EQ.1)GOTO 377
  594. IF(DMAX.LE.xzprec*dmaxge) DMAX = DMAXGE
  595. * facteur de normalisation cf PV pour ne pas avoir de pivot nul
  596. DMAX=DMAX*1.5
  597. * on penalise la matrice en cas de resolution iterative
  598. **pv if (nucrou.eq.1) DMAX=DMAX*1D5
  599. ** write (6,*) ' assem2 i iel dmax ',i,iel,dmax
  600. * XMATRI=IMATTT(IEL)
  601. * SEGACT,XMATRI
  602. c ... dmaxy = maximum des termes dans la première colonne (les 2 premiers exclus) ...
  603. DMAXY=SQRT(XPETIT)*1D5
  604. if (n1.lt.3) then
  605. do ico=1,n1
  606. DMAXY = MAX ( DMAXY, ABS( RE(ICO,1,iel)))
  607. enddo
  608. endif
  609. DO 821 ICO=3,N1
  610. DMAXY = MAX ( DMAXY, ABS( RE(ICO,1,iel)))
  611. 821 CONTINUE
  612. c
  613. IF( IIMPI. EQ.9022 ) WRITE(IOIMP,7398) DMAX
  614. 7398 FORMAT('FACTEUR MULTIPLICATIF DE NORME=',E12.5)
  615. c
  616. c ... copie de la matrice élémentaire fois dmax*coer dans ra ...
  617. DO 21 ICO=1,N1
  618. DO 2110 IKO=1,N1
  619. RA(ICO,IKO)=RE(ICO,IKO,iel)*DMAX*COER
  620. 2110 CONTINUE
  621. 21 CONTINUE
  622. c ... la sous-matrice de dimension 2 est multipliée par dmaxy ...
  623. * commenté car pose un problème dans excite et ne sert sans doute à rien
  624. * write (6,*) ' dmaxy dans assem5 ', dmaxy
  625. xcoef= 1D0
  626. RA(1,1)=RA(1,1)*xcoef
  627. RA(2,1)=RA(2,1)*xcoef
  628. RA(1,2)=RA(1,2)*xcoef
  629. RA(2,2)=RA(2,2)*xcoef
  630. c ... mise à dmax des dnor correspondant aux multiplicateurs ...
  631. DO 22 ICO=1,2
  632. c-old dnor(ival(ico))=dmax
  633. DNOR(IVAL(ICO))=DNOR(IVAL(ICO))*DMAX
  634. 22 CONTINUE
  635. c
  636. c ... boucle sur les variables primales ...
  637. DO 24 ICO=1,N3
  638. c ... ino = n° local du noeud ...
  639. INO=INUINV(NUM(NOELEP(ICO),IEL))
  640. c ... io = n° de ligne du ddl n° ico ...
  641. IO=IVAL(ICO)
  642. if (ico.eq.1) io1=io
  643. if (ico.eq.2) io2=io
  644. c ... cas des noeuds non totalement maitres ...
  645. IF(IO.GT.NBNNMA) THEN
  646. DO 1161 JJ=1,NTTMAI(/1)
  647. IF (NTTMAI(JJ).EQ.INO) GOTO 1171
  648. 1161 CONTINUE
  649. GOTO 1181
  650. 1171 CONTINUE
  651. INO=NNOE+JJ
  652. 1181 CONTINUE
  653. ENDIF
  654. c ... llign contient les lignes liées aux noeud ino ...
  655. LLIGN=ILIGN(INO)
  656. SEGACT,LLIGN*MOD
  657. c ... on rajoute le terme diagonal au diag ...
  658. DIAG(IO)=DIAG(IO)+RA(ICO,ICO)
  659. c ... immmm contient les n°s des lignes des ddl ...
  660. DO 132 JLIJ=1,IMMMM(/1)
  661. JLIJ1=JLIJ
  662. c ... on est censé trouver le bon n° de ligne et aller au 133 ...
  663. IF( IMMMM(JLIJ).EQ.IO) GO TO 133
  664. 132 CONTINUE
  665. c
  666. IF(IIMPI.EQ.1524) WRITE(IOIMP,7354)
  667. 7354 FORMAT( ' PREMIERE ERREUR 5')
  668. CALL ERREUR(5)
  669. RETURN
  670. c
  671. 133 CONTINUE
  672. c
  673. c ... deuxième niveau de boucle sur les variables primales ...
  674. DO 26 IRO=1,N3
  675. c ... ia = n° de colonne du ddl n° iro ...
  676. IA=IVAL(IRO)
  677. c ... cas symétrique !!! ...
  678. IF(IA.GT.IO) GO TO 26
  679. c
  680. c ... jld et jlt sont égaux au début et à la fin (dans xxva et linc)
  681. c des données relatives au ddl n° iro ...
  682. JLT=IPPO(JLIJ1+1)
  683. JLD=IPPO(JLIJ1)+1
  684. c
  685. c ... on cherche le bon n° de colonne dans linc ...
  686. DO 134 JL=JLD,JLT
  687. JL1=JL
  688. IF(LINC(JL).EQ.IA) GO TO 135
  689. 134 CONTINUE
  690. c
  691. IF(IIMPI.EQ.1524) WRITE(IOIMP,7355)
  692. 7355 FORMAT( ' DEUXIEME ERREUR 5')
  693. CALL ERREUR(5)
  694. RETURN
  695. c
  696. 135 CONTINUE
  697. c ... on rajoute le coefficient dans ra au xxva correspondant ...
  698. XXVA(JL1)=XXVA(JL1)+RA(ICO,IRO)
  699. 26 CONTINUE
  700. SEGDES,LLIGN
  701. 24 CONTINUE
  702. * on stocke dans ittr les couples de LX
  703. ittr(io1)=io2
  704. ittr(io2)=io1
  705. * SEGDES,XMATRI
  706. 14 CONTINUE
  707. c ... après son assemblage, une matrice de blocage n'est plus stigmatisée comme telle ...
  708. INOMUL(I)=.TRUE.
  709. 377 CONTINUE
  710. SEGSUP,RA
  711. 11 CONTINUE
  712. c ... fin de la boucle sur les rigidités ...
  713. GO TO 375
  714. c ... fin de la boucle cachée ...
  715. 3750 CONTINUE
  716. c
  717. DO 18 IK=1,NNR
  718. INCTRA=INCTRR(IK)
  719. SEGSUP,INCTRA
  720. 18 CONTINUE
  721. * pv on desactive tout
  722. NNR=IRIGEL(/2)
  723. DO 2 IRI=1,NNR
  724. DESCR=IRIGEL(3,IRI)
  725. SEGDES DESCR
  726. IPT1=IRIGEL(1,IRI)
  727. SEGDES IPT1
  728. xMATRI=IRIGEL(4,IRI)
  729. SEGDES xMATRI
  730. 2 CONTINUE
  731. INTERR(1)=NJTOT
  732. c
  733. IF(IIMPI.GE.1) WRITE(IOIMP,4821) LLVNUL,NJTOT
  734. 4821 FORMAT('Nbre de valeurs non nulles à l assemblage ',I9,/
  735. & 'Nbre de valeurs à l assemblage ',I9)
  736. c
  737. c
  738. SEGDES,SNTT
  739. SEGSUP,INCTRR
  740. SEGDES,MDIAG
  741. SEGDES,MIMIK
  742. SEGDES,MDNOR
  743. SEGSUP,ITOPO
  744. SEGSUP,IITOP
  745. SEGDES,MILIGN
  746. SEGSUP,INUINV
  747. SEGDES,MMATRI
  748. MMMTRI=MMATRI
  749. SEGDES,MINCPO
  750. SEGSUP,IVAL,ITRA,JNOMUL,vmax
  751. IF(NORINC.NE.0) THEN
  752. CALL ASSEM0(MRIGID,2,INWUIT)
  753. SEGDES,MRIGID
  754. ELSE
  755. SEGDES,MRIGID
  756. ENDIF
  757. RETURN
  758. END
  759. c
  760. c
  761. c
  762. c
  763. c
  764.  
  765.  
  766.  
  767.  
  768.  
  769.  
  770.  
  771.  
  772.  
  773.  
  774.  
  775.  
  776.  
  777.  
  778.  
  779.  
  780.  
  781.  
  782.  
  783.  
  784.  
  785.  
  786.  
  787.  
  788.  
  789.  

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