Télécharger assem5.eso

Retour à la liste

Numérotation des lignes :

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