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  24/10/22    21:15:02     12043          
  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.       REAL*8 DMAX,COER,DDDD,DMAXY
  71.       LOGICAL NOMUL
  72. *
  73. * en cas de normalisation des variables.
  74. *
  75.       INWUIT=0
  76.       IF(NORINC.NE.0) CALL ASSEM0(ITRAV1,1,INWUIT)
  77. *
  78. *  pv  on active une fois pour toutes les meleme descr... de la rigidite
  79. *      on en profite pour creer inomul
  80. c
  81. c  ****  recherche de la dimension max de ival,et segini de ival et itra
  82. c
  83. c
  84.       INCTRR=INCTR1
  85.       SEGACT,INCTRR
  86. c
  87.       MRIGID=ITRAV1
  88.  
  89.       SEGACT,MRIGID
  90.       NNR=IRIGEL(/2)
  91.       NNN=0
  92.       SEGINI JNOMUL
  93.       DO 1 IRI=1,NNR
  94. c
  95.           INCTRA=INCTRR(IRI)
  96.           SEGACT INCTRA
  97. c
  98.           DESCR=IRIGEL(3,IRI)
  99.           SEGACT, DESCR
  100. c
  101.           IPT1=IRIGEL(1,IRI)
  102.           SEGACT IPT1
  103. c
  104.           XMATRI=IRIGEL(4,IRI)
  105.           SEGACT XMATRI
  106. c
  107. c    ... na = nb de variables primales ...
  108.           NA=LISINC(/2)
  109. c    ... nnn = max des nb de variables primales ...
  110.           NNN=MAX(NA,NNN)
  111. c
  112. c    ... inomul (non multiplicateurs ?) dit s'il ne s'agit pas des cl ...
  113.           INOMUL(IRI)=.TRUE.
  114.           IF(IPT1.ITYPEL.EQ.49) INOMUL(IRI)=.FALSE.
  115. c
  116.    1  CONTINUE
  117. c
  118. c ... ival a la taille nnn ...
  119.       SEGINI,IVAL
  120. c ... itra a la taille (nnn,2) ...
  121.       SEGINI,ITRA
  122. c
  123. c  ****  activation des segments de travails et de mmatri
  124. c
  125. c ... recherche du nombre d'inconnues impliquees
  126. c
  127.       MMATRI=MMMTRI
  128.       SEGACT,MMATRI*MOD
  129. c
  130.       MINCPO=IINCPO
  131.       SEGACT,MINCPO
  132. c
  133. c
  134.       MAXIPO=0
  135.       DO  I=1,INCPO(/2)
  136.           DO  J=1,INCPO(/1)
  137.               MAXIPO=MAX(INCPO(J,I),MAXIPO)
  138.           ENDDO
  139.       ENDDO
  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.       ntsup=0  
  235. c
  236. c ... la taille de mdnor est inc (nombre d'équations) ...
  237.       SEGINI MDNOR
  238.       IDNORM=MDNOR
  239.       DO 20 I=1,INC
  240.           DNOR(I)=1.D0
  241.    20 CONTINUE
  242. c
  243. c ... si une normalisation est imposée, les valeurs dans mdnor seront modifiées
  244. c    par assem3 ...
  245.       IF(NORINC.NE.0) CALL ASSEM3(MDNOR,MIMIK,MINCPO)
  246. c
  247. c  ****  boucle   *100* sur les numeros de noeuds que l'on assemble
  248. c        -------------------------------------------------------------
  249. c
  250. c    ... ipre = n° de la première ligne concernée par le noeud ino ...
  251.       IPRE=1
  252. c
  253.       DO 999 IJK=1,2
  254. c
  255.           DO 100 INO=1,NNOE
  256.               DO 101 IIT=1,INCDIF
  257.                   MTRA(IIT)=0
  258.  101          CONTINUE
  259. c    ... ider = n° de la dernière ligne concernée par le noeud ino ...
  260.               IDER=0
  261. c
  262.               IF(IJK.EQ.2) IPRE=MAXIPO+1
  263.               DO 998 II=1,INCDIF
  264.                   IF (IJK.EQ.1.AND.INCPO(II,INO).LE.NBNNMA)
  265.      #            IDER=MAX(INCPO(II,INO),IDER)
  266.                   IF (IJK.EQ.2.AND.INCPO(II,INO).GT.NBNNMA) THEN
  267.                       IDER=MAX(INCPO(II,INO),IDER)
  268.                       IPRE=MIN(INCPO(II,INO),IPRE)
  269.                   ENDIF
  270.  998          CONTINUE
  271.               IF (IDER.EQ.0.OR.IPRE.EQ.(MAXIPO+1)) GOTO 100
  272.               LLVVA=0
  273. c
  274. c  ****  boucle   *99*   sur les elements touchant le noeud ino
  275. c                        pour les elements multiplicateur on ne fait pas
  276. c                        l'assemblage
  277. c
  278. c    ... maxele = nb d'éléments qui contiennent le noeud no ino ...
  279.               MAXELE= (IITOP(INO+1) -IITOP(INO))/2
  280. c
  281.  
  282.               DO 99 IELE=1,MAXELE
  283.                   IIU=IITOP(INO) + IELE + IELE - 2
  284. c       ... iel = numero de l'élément dans son maillage ...
  285.                   IEL=ITOPO(IIU)
  286. c       ... iri = numero du maillage (dans irigel) de cet élément ...
  287.                   IRI=ITOPO(IIU+1)
  288. c       ... meleme = pointeur vers ce maillage ...
  289.                   MELEME=IRIGEL(1,IRI)
  290. c       ... descr = pointeur vers le segment descr concerné ...
  291.                   DESCR=IRIGEL(3,IRI)
  292. c       ... inctra contient les indices (dans imik et ihar) des descriptions
  293. c           des ddl correspondants ...
  294.                   INCTRA=INCTRR(IRI)
  295. c       ... xmatri = tableau de pointeurs vers les matrices élémentaires ...
  296.                   XMATRI=IRIGEL(4,IRI)
  297.  
  298. c       ... coer = coefficient multiplicateur ...
  299.                   COER=COERIG(IRI)
  300. c
  301. c  ****  nomul =.false.  il existe  un multiplicateur
  302. c  ****  initialisation de ival. ival(i)=j veut dire que
  303. c  ****  la i eme ligne de la petite matrice s'assemble dans
  304. c  ****  la j eme de la grande matrice.
  305. c
  306. c       ... nin = nombre de variables primales (dans descr) ...
  307.                   NIN=LISINC(/2)
  308. c       ... non multiplicateur ? ...
  309.                   NOMUL=INOMUL(IRI)
  310.                   NA=0
  311. c       ... boucle sur les variables primales de l'élément ...
  312. c
  313.                   DO 98 ICO=1,NIN
  314. c          ... ija = numéro local du noeud-support du ddl ...
  315.                       IJA=INUINV(NUM(NOELEP(ICO),IEL))
  316. c          ... ijb = indice du ddl ...
  317.                       IJB=INCTRA(ICO)
  318. c          ... on met dans ival(ico) le n° de l'équation correspondant au noeud ija
  319. c              et ddl no ijb ...
  320.                       IVAL(ICO)=INCPO(IJB,IJA)
  321. c          ... si on a trouvé le bon noeud ...
  322.                       IF(IJA.EQ.INO) THEN
  323. c             ... on incrémente na (le nombre de ddl connectés au noeud ino ?) ...
  324.                           NA=NA+1
  325.                           ITRA(NA,1)=ICO
  326.                           ITRA(NA,2)=IVAL(ICO)
  327.                       ENDIF
  328.   98              CONTINUE
  329. *                  XMATRI=IMATTT(IEL)
  330. *                  SEGACT,XMATRI
  331. c
  332. c  ****  boucle  *95* sur les inconnues de la petite matrice
  333. c
  334.                               DO 90 INCC=1,NA
  335. c             ... inco = le n° d'équation (ou de ligne) du ddl no incc du noeud ino ...
  336.                           INCO=ITRA(INCC,2)
  337.                           IF (IJK.EQ.1) then
  338.                               IF (INCO.GT.IDER) GOTO 90
  339.                           ENDIF
  340. c             ... iloc = n° de ligne relatif / ipre du noeud ino ...
  341. c
  342.                           ILOC=INCO-IPRE+1
  343.                           IF (IJK.EQ.2) THEN
  344.                               IF (ILOC.LE.0) GOTO 90
  345.                           ENDIF
  346. c
  347. c             ... jj = n° (dans l'élément) de la colonne (et de la ligne) correspondant au
  348. c                 ddl n° incc du noeud ino ...
  349.                           JJ=ITRA(INCC,1)
  350.                   DO 95 IK=1,NIN
  351. c          ... io = n° de l'équation correspondant au ddl no ik de l'élément (n° de la colonne) ...
  352.                       IO=IVAL(IK)
  353. c          ... s'il dépasse le dernier intéressant on ne fait rien ...
  354.                       IF(IO.GT.IDER) GO TO 95
  355. c          ... boucle sur les ddl's du noeud ino ...
  356. c
  357. c
  358. c             ... si io > inco on ne fait plus rien dans cette boucle ...
  359. c             ... car c'est le cas symétrique !!! ...
  360.                           IF (IO.GT.INCO)  GO TO 95
  361. c
  362. c             ... ci-dessous on stockera les quelques lignes de la rigidité concernant le
  363. c                 noeud ino dans le segment tratra. les coefficients vont dans xtra, le deuxième
  364. c                 indice (iloc) est le numéro relatif de la ligne, le premier (imtt) n'a rien à voir avec le
  365. c                 numéro de la colonne. celui-ci est stocké dans ntra(imtt,iloc). l'information inverse
  366. c                 (imtt en fonction de io) se trouve dans ltra(io,iloc). on met dans mtra le nombre
  367. c                 de termes différents assemblés par ligne ...
  368.                           ILTT= LTRA(IO,ILOC)
  369.                           IF(ILTT.EQ.0) THEN
  370. c                ... llvva est initialisé à 0 au début de la boucle 100 ...
  371.                               LLVVA=LLVVA+1
  372. c                ... les mtra sont mis à 0 au début de la boucle 100 ...
  373.                               IMMTT=MTRA(ILOC)+1
  374.                               MTRA(ILOC)=IMMTT
  375.                               XTRA(IMMTT,ILOC)=0.D0
  376.                               NTRA(IMMTT,ILOC)=IO
  377.                               LTRA(IO,ILOC)=IMMTT
  378.                               ILTT=IMMTT
  379.                           ENDIF
  380. c
  381. c             ... si ce n'est pas une condition aux limites, on assemble ...
  382.                           IF(NOMUL) THEN
  383. cfaux !!!            xtra(iltt,iloc)=xtra(iltt,iloc)+re(ik,jj)*coer
  384. c                ... cette ligne a été remplacée par ce qui suit car pour les termes
  385. c                    de couplage entre les ddl du même noeud, les termes au-dessus de
  386. c                    la diagonale étaient pris. Ça marchait car les matrices étaient
  387. c                    symétriques, comme maintenant on travaillera aussi sur des asymétriques ...
  388. ***                           IF(JJ.LT.IK) THEN
  389. ***                               IIILIG=IK
  390. ***                               IIICOL=JJ
  391. ***                           ELSE
  392. ***                               IIILIG=JJ
  393. ***                               IIICOL=IK
  394. ***                           ENDIF
  395. c                ... attention ! dans xtra 1er indice est lié à la colonne, deuxième est le numéro
  396. c                    de la ligne, tandis que dans re c'est l'inverse ...
  397.  
  398. *** on exploite la symetrie de re - XTRA comprend COER
  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.   126         CONTINUE
  482.       DO 127 JPA=IPRE,IDER
  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.       if (ijk.eq.2) ntsup = ntsup + jpa - max(nbnnma,jpa-lpa)
  490.  127  continue
  491. c
  492. c    ... njtot = profil total ...
  493.               NJTOT=NJTOT+NJMAX
  494. c    ... ijmax = (profil / noeud) maxi ...
  495.               IF(IJMAX.LT.NJMAX) IJMAX=NJMAX
  496.               SEGDES,LLIGN
  497.               IPRE=IDER+1
  498.   100     CONTINUE
  499.       if (ijk.eq.1) ntint=njtot                   
  500.       iok=1
  501.       if (ntsup.gt.ntint) iok=0
  502. c
  503.  999  CONTINUE
  504. **    write(6,*) 'assem5 surface eliminee surface gardee',ntint,ntsup
  505. **    write(6,*) 'super element ok',iok
  506. c
  507. c ... fin de la boucle sur les noeuds ...
  508.       SEGSUP TRATRA
  509. C
  510. C  ****  ON REPREND TOUTE LES MATRICES CONTENANT LES MULTIPLICATEURS
  511. C  ****  POUR MULTIPLIER TOUS LEURS TERMES PAR UNE NORME ATTACHEE
  512. C  ****  A CHAQUE MULTIPLICATEUR.PUIS ON LES ASSEMBLE.
  513. C
  514. *  d'abord etablir une norme generale pour le cas ou on n'arrive pas
  515. *  a calculer la norme particuliere
  516.       DMAXGE=xpetit
  517.       DO 378 I=1,INC
  518. **    write (6,*) ' assem5 diag vmav ',diag(i),vmax(i)
  519.       DMAXGE=MAX(DMAXGE,abs(vmax(i)))
  520.   378 CONTINUE
  521.       if (dmaxge.lt.xpetit/xzprec*10) dmaxge=1.d0
  522.       if (iimpi.eq.1524)
  523.      >  write (6,*) ' nb inconnues facteur multiplicatif general ',
  524.      >  INC,DMAXGE
  525.       IENMU=0
  526. c ... attention ! ici commence une boucle cachée (exécutée avec un goto 375) ...
  527.   375 IENMU1 = IENMU
  528.       IENMU=0
  529. c ... boucle sur les rigidités qui calcule le nombre de matrices de blocages ...
  530.       DO 376 I=1,NNR
  531.           IF(.NOT.INOMUL(I)) IENMU=IENMU+1
  532.  376  CONTINUE
  533. c
  534. c ... s'il n'y en a pas on saute cette partie du code ...
  535.       IF( IENMU.EQ.0) GO TO 3750
  536. c
  537. c ... mimik contient les noms des variables primales ...
  538.       MIMIK=IIMIK
  539.       SEGACT,MIMIK
  540. c ... boucle sur les rigidités ...
  541.       DO 11 I=1,NNR
  542. c    ... si celle si n'est pas une matrice de blocage on passe à la suivante ...
  543.           IF(INOMUL(I)) GO TO 11
  544. c
  545.           DESCR=IRIGEL(3,I)
  546. c    ... n3 = nb de variables primales ...
  547.           N3=LISINC(/2)
  548.           COER=COERIG(I)
  549.           MELEME=IRIGEL(1,I)
  550.           INCTRA=INCTRR(I)
  551.           XMATRI=IRIGEL(4,I)
  552. c    ... n2 = nombre d'éléments dans le support géométrique ...
  553.           N2=NUM(/2)
  554. c    ... À quoi correspond ce cas ? (pas de matrices élémentaires) ...
  555.           IF (RE(/3).EQ.0) THEN
  556.               INOMUL(I)=.TRUE.
  557.               GOTO 11
  558.           ENDIF
  559. *          XMATRI=IMATTT(1)
  560. *          SEGACT,XMATRI
  561. c    ... n1 = nombre de variables duales ...
  562. c    ... pourquoi va-t-on chercher n3 dans descr et n1 dans re ? ...
  563.           N1=RE(/1)
  564.           SEGINI,RA
  565. c    ... boucle sur les éléments ...
  566.           DO 14 IEL=1,N2
  567. c       ... boucle sur les inconnues ...
  568.               DO 15 ICO=1,N3
  569. c          ... ija = n° local du noeud-support de l'inconnue n° ico ...
  570.                   IJA=INUINV(NUM(NOELEP(ICO),IEL))
  571. c          ... ijb = n° de l'inconnue ...
  572.                   IJB=INCTRA(ICO)
  573. c          ... ival contient la correspondance entre le n° local du ddl
  574. c              et le n° d'équation (de ligne) correspondant ...
  575.                   IVAL(ICO)=INCPO(IJB,IJA)
  576.    15         CONTINUE
  577.  
  578.                DDDD= xpetit
  579. c       ... dmax = max des valeurs absolues des termes diagonaux correspondants
  580. c           au ddl de l'élément n° iel ...
  581. c           COER a deja ete applique sur vmax
  582.               DMAX=DDDD
  583. *  la boucle suivante demarre 3 car les deux premiers sont les multiplicateurs de lagrange
  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 (norinc.eq.0) dmaxy=1.d0
  605. *  demarrage a 3 aussi. On a toujours 1 -1 sur les LX
  606.               DO 821 ICO=3,N1
  607.                   DMAXY = MAX ( DMAXY, ABS( RE(ICO,1,iel)))
  608. 821           CONTINUE
  609. c
  610.               IF( IIMPI. EQ.9022 ) WRITE(IOIMP,7398) DMAX
  611.  7398         FORMAT('FACTEUR MULTIPLICATIF DE NORME=',E12.5)
  612. ** coer est dans dmax, mais pas dans dmaxy
  613.       DMAX = DMAX / DMAXY       
  614.  
  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)*coer*DMAX  
  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.             if (norinc.eq.0.d0) dmaxy=dmaxy*2.d0
  626.             RA(1,1)=RA(1,1)*DMAXY
  627.             RA(2,1)=RA(2,1)*DMAXY
  628.             RA(1,2)=RA(1,2)*DMAXY
  629.             RA(2,2)=RA(2,2)*DMAXY
  630. c       ... mise à dmax des dnor correspondant aux multiplicateurs ...
  631.               DO 22 ICO=1,2
  632.                   dnor(ival(ico))=dmax   
  633.    22         CONTINUE
  634. c
  635. c       ... boucle sur les variables primales ...
  636.               DO 24 ICO=1,N3
  637. c          ... ino = n° local du noeud ...
  638.                   INO=INUINV(NUM(NOELEP(ICO),IEL))
  639. c          ... io = n° de ligne du ddl n° ico ...
  640.                   IO=IVAL(ICO)
  641.       if (ico.eq.1) io1=io
  642.       if (ico.eq.2) io2=io
  643. c          ... cas des noeuds non totalement maitres ...
  644.                   IF(IO.GT.NBNNMA) THEN
  645.                       DO 1161 JJ=1,NTTMAI(/1)
  646.                           IF (NTTMAI(JJ).EQ.INO) GOTO 1171
  647.  1161                 CONTINUE
  648.                       GOTO 1181
  649.  1171                 CONTINUE
  650.                       INO=NNOE+JJ
  651.  1181                 CONTINUE
  652.                   ENDIF
  653. c          ... llign contient les lignes liées aux noeud ino ...
  654.                   LLIGN=ILIGN(INO)
  655.                   SEGACT,LLIGN*MOD
  656. c          ... on rajoute le terme diagonal au diag ...
  657.                   DIAG(IO)=DIAG(IO)+RA(ICO,ICO)
  658. c          ... immmm contient les n°s des lignes des ddl ...
  659.                   DO 132 JLIJ=1,IMMMM(/1)
  660.                       JLIJ1=JLIJ
  661. c             ... on est censé trouver le bon n° de ligne et aller au 133 ...
  662.                       IF( IMMMM(JLIJ).EQ.IO) GO TO 133
  663.   132             CONTINUE
  664. c
  665.                   IF(IIMPI.EQ.1524) WRITE(IOIMP,7354)
  666.  7354             FORMAT( ' PREMIERE ERREUR 5')
  667.                   CALL ERREUR(5)
  668.                   RETURN
  669. c
  670.   133             CONTINUE
  671. c
  672. c          ... deuxième niveau de boucle sur les variables primales ...
  673.                   DO 26 IRO=1,N3
  674. c             ... ia = n° de colonne du ddl n° iro ...
  675.                       IA=IVAL(IRO)
  676. c             ... cas symétrique !!! ...
  677.                       IF(IA.GT.IO) GO TO 26
  678. c
  679. c             ... jld et jlt sont égaux au début et à la fin (dans xxva et linc)
  680. c                 des données relatives au ddl n° iro ...
  681.                       JLT=IPPO(JLIJ1+1)
  682.                       JLD=IPPO(JLIJ1)+1
  683. c
  684. c             ... on cherche le bon n° de colonne dans linc ...
  685.                       DO 134 JL=JLD,JLT
  686.                           JL1=JL
  687.                           IF(LINC(JL).EQ.IA) GO TO 135
  688.  134                  CONTINUE
  689. c
  690.                       IF(IIMPI.EQ.1524) WRITE(IOIMP,7355)
  691.  7355                 FORMAT( ' DEUXIEME ERREUR 5')
  692.                       CALL ERREUR(5)
  693.                       RETURN
  694. c
  695.  135                  CONTINUE
  696. c             ... on rajoute le coefficient dans ra au xxva correspondant ...
  697.                       XXVA(JL1)=XXVA(JL1)+RA(ICO,IRO)
  698.    26             CONTINUE
  699.                   SEGDES,LLIGN
  700.    24         CONTINUE
  701. *  on stocke dans ittr les couples de LX
  702.       ittr(io1)=io2
  703.       ittr(io2)=io1
  704. *              SEGDES,XMATRI
  705.    14     CONTINUE
  706. c    ... après son assemblage, une matrice de blocage n'est plus stigmatisée comme telle ...
  707.           INOMUL(I)=.TRUE.
  708.   377     CONTINUE
  709.           SEGSUP,RA
  710.    11 CONTINUE
  711. c ... fin de la boucle sur les rigidités ...
  712.       GO TO 375
  713. c ... fin de la boucle cachée ...
  714.  3750 CONTINUE
  715. c
  716.       DO 18 IK=1,NNR
  717.           INCTRA=INCTRR(IK)
  718.           SEGSUP,INCTRA
  719.    18 CONTINUE
  720. *  pv on desactive tout
  721.       NNR=IRIGEL(/2)
  722.       DO 2 IRI=1,NNR
  723.           DESCR=IRIGEL(3,IRI)
  724.           SEGDES DESCR
  725.           IPT1=IRIGEL(1,IRI)
  726.           SEGDES IPT1
  727.           xMATRI=IRIGEL(4,IRI)
  728.           SEGDES xMATRI
  729.    2  CONTINUE
  730.       INTERR(1)=NJTOT
  731. c
  732.       IF(IIMPI.GE.1) WRITE(IOIMP,4821) LLVNUL,NJTOT
  733.  4821 FORMAT('Nbre de valeurs non nulles à l assemblage      ',I9,/
  734.      &       'Nbre de valeurs à l assemblage                 ',I9)
  735. c
  736. c
  737.       SEGDES,SNTT
  738.       SEGSUP,INCTRR
  739.       SEGDES,MDIAG
  740.       SEGDES,MIMIK
  741.       SEGDES,MDNOR
  742.       SEGSUP,ITOPO
  743.       SEGSUP,IITOP
  744.       SEGDES,MILIGN
  745.       SEGSUP,INUINV
  746.       SEGDES,MMATRI
  747.       MMMTRI=MMATRI
  748.       SEGDES,MINCPO
  749.       SEGSUP,IVAL,ITRA,JNOMUL,vmax
  750.       IF(NORINC.NE.0) THEN
  751.           CALL ASSEM0(MRIGID,2,INWUIT)
  752.           SEGDES,MRIGID
  753.       ELSE
  754.           SEGDES,MRIGID
  755.       ENDIF
  756.       RETURN
  757.       END
  758. c
  759. c
  760. c
  761. c
  762. c
  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.  
  790.  
  791.  
  792.  
  793.  
  794.  
  795.  

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