Télécharger la plaquette Cast3M


fr en

Télécharger assem5.eso

Retour à la liste

Numérotation des lignes : 

assem5

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

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