Télécharger la plaquette Cast3M


fr en

Télécharger asns1.eso

Retour à la liste

Numérotation des lignes : 

asns1

  1. C ASNS1     SOURCE    PV        21/09/20    21:15:00     11104          
  2.       SUBROUTINE ASNS1 ( IPOIRI,MMATRX,INUINY,ITOPOY,IMINIY,IPOY,
  3.      & INCTRY,INCTRZ,IITOPY,ITOPOD,IITOPD,IPODD)
  4.  
  5. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  6. C  CE SUBROUTINE SERT A L'ASSEMBLAGE DE MATRICES SYMETRIQUES
  7. C  EN VUE D'UNE INVERSION PAR UNE METHODE DE CROUT
  8. C
  9. C  EN ENTREE:
  10. C  ****  IPOIRI: POINTEUR SUR OBJET MRIGIDITE,NON MODIFIE
  11. C  EN SORTIE:
  12. C  ****  INUINV IMINI ITOPO IPOY INCTRY  SONT DES POINTEURS DES SEGMENTS
  13. C        DE TRAVAIL SERVANT A L'ASSEMBLAGE ILS SONT DETRUITS EN FIN
  14. C        D'ASSEMBLAGE OU DE TRIANGULARISATION
  15. C  ****  MMATRI EST LE POINTEUR DE L'OBJET FUTUR MATRICE TRIANGULARISEE.
  16. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  17.  
  18.       IMPLICIT INTEGER(I-N)
  19.       IMPLICIT REAL*8(A-H,O-Z)
  20. -INC PPARAM
  21. -INC CCOPTIO
  22. -INC SMELEME
  23. -INC SMCOORD
  24. -INC CCREEL
  25.  
  26.       SEGMENT,IMIN(NNOE)
  27.       SEGMENT,IMINB(NNOE)
  28.       SEGMENT ICPR(nbpts)
  29. C
  30. -INC SMRIGID
  31. -INC SMMATRI
  32. C
  33.       SEGMENT,INUINV(NNGLOB)
  34.       SEGMENT,ITOPO(IENNO)
  35.       SEGMENT,IITOP(NNOE+1)
  36.       SEGMENT,ITOPOB(IENNO)
  37.       SEGMENT,IITOPB(NNOE+1)
  38.       SEGMENT,IMINI(INC)
  39.       SEGMENT,IPOS(NNOE1)
  40.       SEGMENT,IPOD(NNOE1)
  41.       SEGMENT,INCTRR(NIRI)
  42.       SEGMENT,INCTRD(NIRI)
  43.       SEGMENT,INCTRA(NLIGRE)
  44.       SEGMENT DIATMP(maxt,NNOE)
  45.       segment strv
  46.        integer itrv1(maxt)
  47.        integer itrv2(maxt)
  48.        real*8 dtrv1(maxt)
  49.        real*8 dtrv2(maxt)
  50.       endsegment
  51.       segment mondu
  52.        character*4 mondua(nnn)
  53.        integer ipris(nnn),inosel(nnn)
  54.       endsegment
  55.  
  56. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  57. C  ****  CES TABLEAUX SERVENT AU REPERAGE DE LA MATRICE POUR L'ASSEMBLAG
  58. C  ****  IL SERONT TOUS SUPPRIMES EN FIN D'ASSEMBLAGE.
  59. C
  60. C
  61. C  ****   MAXINC= MAXIMUM DE COMPOSANTES CONCERNANT UN NOEUD
  62. C
  63. C
  64. C  ****   IITOP(K)=I LE 1ER ELEMENT TOUCHANT LE NOEUD K SE TROUVE EN
  65. C                    IEME POSITION DANS ITOPO
  66. C  ****   ITOPO(I)=L: LE 1ER ELEMENT TOUCHANT LE K EME NOEUD DE LA
  67. C         ITOPO(I+1)=M  MATRICE EST LE LIEME  DE L'OBJET GEOMETRIE
  68. C                       DEFINI PAR LE POINTEUR M
  69. C  ****   IPOS(I)=J   : LA 1 ERE INCONNUE DU NOEUD I EST EN J+1 EME
  70. C                       POSITION
  71. C  ****  IMINI(I)=J     LA PLUS PETITE INCONNUE QUI EST RELIEE A LA IEME
  72. C                       EST L'INCONNUE J.
  73. C  ****  INUINV(I)=J    J EST LE NOUVEAU NUMERO DU NOEUD I
  74. C
  75. C  ****  INCTRR(NIRI) - NIRI=NRIGEL du IPOIRI (objet MRIGID passé en argument)
  76. C                       pointeurs sur INCTRA
  77. C
  78. C  Variables locales :
  79. C  --------------------
  80. C     * NNVA = NRIGEL (nombre d'objets MRIGID élémentaires) dans IPOIRI (objet
  81. C              MRIGID) passé en argument)
  82. C     * NLIGRE = NLIGRP - nombre de variables primales (dans un segment DESCR)
  83. C     * IMELP = pointeur d'un MELEME contenant un noeud "normal"
  84. C     * ICDOUR = ???
  85. C
  86. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  87.  
  88.       CHARACTER*4 CNOHA,lisi
  89.       integer*4   noha
  90.       equivalence (cnoha,noha)
  91.       DATA CNOHA/'NOHA'/
  92.       DATA IPOIN/1/
  93.  
  94.       NNGLOB=nbpts
  95.  
  96.       MRIGID=IPOIRI
  97.       SEGACT,MRIGID
  98.  
  99.       NNVA=IRIGEL(/2)
  100.       NIRI=NNVA
  101.       SEGINI,INCTRR
  102.       SEGINI,INCTRD
  103.  
  104.       MELEME=IRIGEL(1,1)
  105.       SEGACT MELEME
  106. C ... ITYPEL = 27 correspond aux éléments 'ATTA' ...
  107.       IF(ITYPEL.NE.27) GO TO 801
  108.       SEGDES MELEME
  109. C
  110. C  **** ASSEMBLAGE DANS LE CAS DE L'ANALYSE MODALE. ON COMPTE LES POINTS
  111. C  **** DANS ICPR
  112. C
  113.       SEGINI INUINV,ICPR
  114.       IKI=0
  115.       DO 700 I=1,NNVA
  116.          MELEME=IRIGEL(1,I)
  117.          SEGACT MELEME
  118.          NBNN=NUM(/1)
  119.          NBELEM=NUM(/2)
  120.          DO  I1=1,NBELEM
  121.             DO 701 I2=1,NBNN
  122.                IP1=NUM(I2,I1)
  123.                IF(ICPR(IP1).NE.0) GO TO 701
  124.                IKI=IKI+1
  125.                ICPR(IP1)=IKI
  126.  701     CONTINUE
  127.          enddo
  128.          SEGDES MELEME
  129.  700  CONTINUE
  130. C
  131. C  **** FABRICATION DU TABLEAU INUINV
  132. C  **** ON MET LES POINTS QUI ONT POUR INCONNUE ALFA EN TETE
  133. C
  134.       NNOE=IKI
  135.       NBETA=0
  136.       DO 710 I=1,NNVA
  137.          MELEME=IRIGEL(1,I)
  138.          DESCR =IRIGEL(3,I)
  139.          SEGACT MELEME,DESCR
  140.          NBNN=NUM(/1)
  141.          NBELEM=NUM(/2)
  142.          NLIGRE=LISINC(/2)
  143.          DO  I1=1,NBELEM
  144.             DO 711 I2=1,NBNN
  145.                IP1=NUM(I2,I1)
  146.                IF(ICPR(IP1).EQ.0) GO TO 711
  147.  715           CONTINUE
  148.                NBETA=NBETA+1
  149.                IKI=NNOE-NBETA+1
  150.  716           CONTINUE
  151.                INUINV(IP1)=IKI
  152.                ICPR(IP1)=0
  153.  711     CONTINUE
  154.          enddo
  155.          SEGDES MELEME,DESCR
  156. *        SEGSUP IPB
  157.  710  CONTINUE
  158.       SEGSUP ICPR
  159.       ICDOUR=NNOE
  160.       GO TO 800
  161. C
  162. C  ****  ON FABRIQUE UN NOUVEL OBJET GEOMETRIE CONTENANT TOUTES LES
  163. C  ****  GEOMETRIES ELEMENTAIRES. CET OBJET CONTIENT NNVA OBJETS
  164. C  ****  GEOMETRIQUES ELEMENTAIRES. PUIS ON ENVOIE DANS NUMOPT QUI
  165. C  ****  FOURNIT EN RETOUR INUINV(NUM(I,J))=K DONNE LE NOUVEAU
  166. C  ****  NUMERO LOCAL DU POINT NUM(I,J).K VARIE DE 1 A ICDOUR.
  167. C  ****  LE PREMIER NOEUD DE L'OBJET GEOMETRIQUE EST LE PREMIER NOEUD
  168. C  ****  DE LA MATRICE, ETC...
  169. C
  170.  801  CONTINUE
  171.       IKK=1
  172.  722  CONTINUE
  173.       MELEME=IRIGEL(1,IKK)
  174.       SEGACT,MELEME
  175.  
  176.       DESCR=IRIGEL(3,IKK)
  177.       SEGACT,DESCR
  178.  
  179.       NLIGRE=LISINC(/2)
  180.       DO 720 K=1,NLIGRE
  181.          IF(LISINC(K).NE.'LX  ') GO TO 721
  182.   720 CONTINUE
  183.       SEGDES,MELEME
  184.       SEGDES,DESCR
  185.       IKK=IKK+1
  186.       IF(IKK.LE.NNVA)  GO TO 722
  187.       DO 4862 I=1,NNVA
  188.          MELEME= IRIGEL(1,I)
  189.          SEGACT MELEME
  190.          IF(ITYPEL.EQ.49) THEN
  191.             DESCR=IRIGEL(3,I  )
  192.             SEGACT,DESCR
  193.             K = 3
  194.             GO TO 721
  195.          ELSE
  196.             SEGDES MELEME
  197.          ENDIF
  198.  4862 CONTINUE
  199.       K=1
  200.       MELEME= IRIGEL(1,1)
  201.       DESCR= IRIGEL(3,1)
  202.       SEGACT MELEME,DESCR
  203.  
  204. C ... On arrive ici si :
  205. C     * LISINC(K) != 'LX  ' => K est le premier parmi K tels que LISINC(K) != 'LX  '
  206. C     * ITYPEL d'un des maillages = 49 (élément 'MULT') => K = 3
  207. C     * tous les autres cas => K = 1
  208.  
  209. C ... IA = numéro (dans l'élément) du noeud concerné par le DDL No K ...
  210. C ... I1 = numéro (absolu) du noeud concerné par le DDL No K,
  211. C     Ce noeud sera mis dans un MELEME dont le pointeur est stocké dans IMELP ...
  212.  
  213.   721 IA=NOELEP(K)
  214.       I1=NUM(IA,1)
  215.       SEGDES,MELEME,DESCR
  216.  
  217.       NBSOUS=0
  218.       NBNN=1
  219.       NBREF=0
  220.       NBELEM=1
  221.       SEGINI,MELEME
  222.       ITYPEL=1
  223.       NUM(1,1)=I1
  224.       SEGDES,MELEME
  225.       IMELP=MELEME
  226.  
  227. C ... Le MELEME créé ici est un MELEME composé qui contiendra le MELEME
  228. C     pointé par IMELP et tous les MELEME pointés par IRIGEL(1,*) ...
  229.       NBSOUS=NNVA+1
  230.       NBREF=0
  231.       NBNN=0
  232.       NBELEM=0
  233.       SEGINI,MELEME
  234.       LISOUS(1)=IMELP
  235.       DO 12 I=1,NNVA
  236.          ipt1=irigel(1,i)
  237.          segact ipt1
  238.           if (irigel(7,i).ne.0.and.ipt1.lisref(/1).ne.0.and.
  239.      >          ipt1.itypel.eq.49) then
  240.           write(6,*) 'assemblage condition non symetrique'
  241.           ipt1=ipt1.lisref(1)
  242.          endif
  243.          LISOUS(I+1)=ipt1
  244.    12 CONTINUE
  245.       SEGDES,MELEME
  246.  
  247.       ICDOUR=0
  248.       SEGINI,INUINV
  249.       SEGDES,INUINV
  250.       CALL NUMOPT(MELEME,INUINV,ICDOUR)
  251. C ... A la sortie INUINV contient l'ordre des noeuds et ICDOUR le nombre de noeuds présents dans MELEME ...
  252.       SEGACT INUINV
  253.       SEGSUP,MELEME
  254.       MELEME=IMELP
  255.       SEGDES,MELEME
  256. C
  257. C  ****  CREATION D'UN OBJET GEOMETRIE QU'IL FAUDRA CHANGER EN CAS DE
  258. C  ****  RENUMEROTATION GENERALE.ON PROFITE DE LA BOUCLE POUR CREE LE
  259. C  ****  TABLEAU IMIN(I)=J QUI DIT QUE J ELEMENTS TOUCHE LE NOEUD I(NU-
  260. C  ****  MEROTATION LOCALE).
  261. C
  262.  800  CONTINUE
  263.       NNOE=ICDOUR
  264.       SEGINI,IMIN,IMINB
  265.       NNOE1=NNOE+1
  266.       SEGINI,IPOS,IPOD
  267.       NBSOUS=0
  268.       NBREF=0
  269.       NBNN=1
  270.       NBELEM=ICDOUR
  271.       SEGINI,IPT1
  272.       IPT1.ITYPEL=IPOIN
  273.       DO 16 IRI=1,NNVA
  274.          DO  I=1,NNOE
  275.            ipod(I)=0
  276.             IPOS(I)=0
  277.          enddo
  278.          MELEME=IRIGEL(1,IRI)
  279.          SEGACT,MELEME
  280.          DESCR=IRIGEL(3,IRI)
  281.          segact descr
  282.          N1=NUM(/1)
  283.          N2=NUM(/2)
  284. *         write(6,*) 'noelep', ( noelep(iu),iu=1,noelep(/1))
  285. *         write(6,*) 'noeled', ( noeled(iu),iu=1,noeled(/1))
  286.          DO 17 I=1,N2
  287.             DO 171 J=1,NOELEP(/1)
  288.              K = NUM( NOELEP(J),I)
  289.              M=INUINV(K)
  290.              IF(IPOS(M).NE.I) THEN
  291.                 IMIN(M)=IMIN(M)+1
  292.                 IPT1.NUM(1,M)=K
  293.                 IPOS(M)=I
  294.              ENDIF
  295.   171      CONTINUE
  296.           DO 172 J=1,NOELED(/1)
  297.              K = NUM( NOELED(J),I)
  298.              M=INUINV(K)
  299.              IF(IPOD(M).NE.I) THEN
  300.                 IMINB(M)=IMINB(M)+1
  301.                 IPOD(M)=I
  302.              ENDIF
  303.   172      CONTINUE
  304.    17    CONTINUE
  305.          SEGDES,MELEME
  306.    16 CONTINUE
  307. C
  308. C  ****  INITIALISATION DE ITOPO. ON UTILISE IMIN POUR SE POSITIONNER
  309. C  ****  DANS ITOPO .
  310. C  ... ITOPO contiendra pour chaque noeud et chaque élément contenant
  311. C      ce noeud 2 nombres :
  312. C            1. numéro de l'élément dans son maillage
  313. C            2. numéro du maillage (dans IRIGEL) de cet élément
  314. C
  315. C ... IITOP servira pour déterminer la taille de ITOPO ainsi que pour
  316. C     se retrouver dedans ...
  317. C
  318.       SEGINI,IITOP,IITOPB
  319.       IITOP(1)=1
  320.       IITOPB(1)=1
  321. *       write(6,*) ' imin', ( imin(iu),iu=1,imin(/1))
  322. *       write(6,*) ' iminb', ( iminb(iu),iu=1,iminb(/1))
  323.       DO 18 I=1,NNOE
  324.          IITOP(I+1)=IMIN(I)* 2  + IITOP(I)
  325.          IITOPB(I+1)=IMINB(I)* 2  + IITOPB(I)
  326.    18 CONTINUE
  327.       DO  I=1,NNOE
  328.          IMINB(I)=0
  329.          IMIN(I)=0
  330.       enddo
  331. C ... IENNO = taille d'ITOPO ...
  332.       IENNO=IITOP(NNOE+1)
  333.       SEGINI,ITOPO
  334.       IENNO=IITOPB(NNOE+1)
  335.       SEGINI ITOPOB
  336.       DO 21 IRI=1,NNVA
  337.          DO I=1,NNOE
  338.             IPOD(I)=0
  339.             IPOS(I)=0
  340.          enddo
  341.          MELEME=IRIGEL(1,IRI)
  342.          SEGACT,MELEME
  343.          DESCR = IRIGEL(3,IRI)
  344.          N2=NUM(/2)
  345.          DO 22 I=1,N2
  346.              DO 221 J=1,NOELEP(/1)
  347.                M=INUINV(NUM(NOELEP(J),I))
  348.                IF(IPOS(M).NE.I) THEN
  349.                   IMIN(M)=IMIN(M)+1
  350.                   IUY= 2* ( IMIN(M)-1 )  + IITOP(M)
  351. C             ... Remplissage d'ITOPO ...
  352.                   ITOPO(IUY)=I
  353.                   ITOPO(IUY+1)=IRI
  354.                   IPOS(M)=I
  355.                ENDIF
  356.   221    CONTINUE
  357.              DO 222 J=1,NOELED(/1)
  358.                M=INUINV(NUM(NOELED(J),I))
  359.                IF(IPOD(M).NE.I) THEN
  360.                   IMINB(M)=IMINB(M)+1
  361.                   IUY= 2* ( IMINB(M)-1 )  + IITOPB(M)
  362. C             ... Remplissage d'ITOPO ...
  363.                   ITOPOB(IUY)=I
  364.                   ITOPOB(IUY+1)=IRI
  365.                   IPOD(M)=I
  366.                ENDIF
  367.   222    CONTINUE
  368.    22 CONTINUE
  369.          SEGDES,MELEME
  370.    21 CONTINUE
  371. C
  372. C  RECHERCHE DE LA VALEUR PAR DEFAUT DE L'HARMONIQUE DANS LE CAS
  373. C  DE L'UTILISATION DE " OPTION MODE FOUR NOHAR "
  374. C
  375. C ... On passe cette boucle sans erreur si tous les IRIGEL(5,*) sont égaux
  376. C     soit à NOHA soit à une autre valeur fixe (IARDEF) ...
  377. C
  378.       DO 230 IRI=1,NNVA
  379.          IHARIR=IRIGEL(5,IRI)
  380.          IF( IHARIR . NE. NOHA) THEN
  381.             IARDEF = IHARIR
  382.             GO TO 231
  383.          ENDIF
  384.   230 CONTINUE
  385.       CALL ERREUR ( 21)
  386.       RETURN
  387.  
  388.   231 CONTINUE
  389.       DO 232 IRI=1,NNVA
  390.          IF( IRIGEL(5,IRI) .EQ.NOHA) GO TO 232
  391.          IF( IRIGEL(5,IRI).EQ.IARDEF ) GO TO 232
  392.          CALL ERREUR (21)
  393.          RETURN
  394.  232  CONTINUE
  395. C
  396. C  ****  RECHERCHE DE LA VALEUR MAXINC QUI PERMET DE DIMENSIONNER INCPOS
  397. C
  398. C ... Les quatre segments sont à l'origine de longueur nulle ...
  399.       SEGINI,MIDUA
  400.       SEGINI,MIMIK
  401.       SEGINI,MHARK
  402.       SEGINI,MHAR1
  403.  
  404.       DESCR=IRIGEL(3,1)
  405.       SEGACT,DESCR
  406.  
  407.       IMIK(**)=LISINC(1)
  408.       IAAR=IRIGEL(5,1)
  409.       IF(IAAR.EQ.NOHA) IAAR = IARDEF
  410.       IHAR(**)= IAAR
  411.       IDUA(**)=LISDUA(1)
  412.       MHAR1.IHAR(**)= IAAR
  413.  
  414.       MAXINC=1
  415.       DO 23 IRI=1,NNVA
  416.          DESCR=IRIGEL(3,IRI)
  417.          IHARIR=IRIGEL(5,IRI)
  418.          IF(IHARIR. EQ.NOHA ) IHARIR = IARDEF
  419.          SEGACT,DESCR
  420.          NLIGRE=LISINC(/2)
  421.          DO 26 I=1,NLIGRE
  422.             DO 24 J=1,MAXINC
  423.                IF(IMIK(J).NE.LISINC(I)) GO TO 24
  424.                IF(IHAR(J).EQ.IHARIR) GO TO 26
  425.    24       CONTINUE
  426. C       ... On empile les valeurs d'IHARIR et LISINC dans
  427. C           leurs segments si le couple (IHARIR,LISINC) n'y est pas
  428. C           encore représenté ...
  429.             MAXINC=MAXINC+1
  430.             IHAR(**)=IHARIR
  431.             IMIK(**)=LISINC(I)
  432.    26    CONTINUE
  433.    23 CONTINUE
  434.  
  435.       MAXDUA=1
  436.       DO 2322 IRI=1,NNVA
  437.          DESCR=IRIGEL(3,IRI)
  438.          IHARIR=IRIGEL(5,IRI)
  439.          IF(IHARIR. EQ.NOHA ) IHARIR = IARDEF
  440.          SEGACT,DESCR
  441.          NLIGRE=LISDUA(/2)
  442.          DO 262 I=1,NLIGRE
  443.             DO 242 J=1,MAXDUA
  444.                IF(IDUA(J).NE.LISDUA(I)) GO TO 242
  445.                IF(MHAR1.IHAR(J).EQ.IHARIR) GO TO 262
  446.   242       CONTINUE
  447. C       ... On empile les valeurs d'IHARIR et LISDUA dans
  448. C           leurs segments si le couple (IHARIR,LISDUA) n'y est pas
  449. C           encore représenté ...
  450.             MAXDUA=MAXDUA+1
  451.             MHAR1.IHAR(**)=IHARIR
  452.             IDUA(**)=LISDUA(I)
  453.   262    CONTINUE
  454.          SEGDES,DESCR
  455.  2322 CONTINUE
  456. *      write(6,*) ' imik'
  457. *      write(6,*) ( imik(iu),iu=1,imik(/2))
  458. *      write(6,*) ' idua avant'
  459. *      write(6,*) ( idua(iu),iu=1,idua(/2))
  460.       nnn = idua(/2)
  461.       nqq = imik(/2)
  462.       if(nnn.ne.nqq) then
  463. *  on verra plus tard
  464.        call erreur(5)
  465.        stop
  466.       endif
  467. * petit travail pour mettre dans le meme ordre les inconnues
  468.       segini mondu
  469.       do 476 iu=1,imik(/2)
  470.       lisi=imik(iu)
  471.       do 477 io=1,idua(/2)
  472.       if( idua(io).eq.lisi) go to 478
  473.   477 continue
  474.       inosel(iu)=1
  475.       go to 476
  476.   478 mondua(iu)= idua(io)
  477.       ipris(io)=1
  478.   476 continue
  479.       do 472 iu=1,inosel(/1)
  480.       if( inosel(iu).eq.0) go to 472
  481.        do 473 io=1,ipris(/1)
  482.         if( ipris(io).eq.1) go to 473
  483.         ipris(io)=1
  484.          mondua(iu)=idua(io)
  485.          go to 472
  486.   473 continue
  487.   472 continue
  488.       do 479 iu=1,idua(/2)
  489.       idua(iu)=mondua(iu)
  490.   479 continue
  491.       segsup mondu
  492. *      write(6,*) ' idua apres'
  493. *      write(6,*) ( idua(iu),iu=1,idua(/2))
  494. C
  495. C  **** INITIALISATION DE INCPOS ET DE INCTRA.
  496. C
  497. C ... Les dimensions des segments MINCPO initialisés ci-dessous sont les
  498. C     suivantes : MAXI = nombre de différentes variables primales (ou duales)
  499. C                 NNOE = nombre de noeuds effectivement présents
  500.       MAXI=MAXINC
  501.       SEGINI,MINCPO
  502.  
  503.       MAXI=MAXDUA
  504.       SEGINI,MIPO1
  505.       maxt=max(maxinc,maxdua)
  506.       SEGINI DIATMP,strv
  507.  
  508.       DO 29 IRI=1,NNVA
  509.          IHARIR=IRIGEL(5,IRI)
  510.          IF(IHARIR.EQ.NOHA ) IHARIR = IARDEF
  511.  
  512.          DESCR=IRIGEL(3,IRI)
  513.          SEGACT,DESCR
  514.  
  515.          NLIGRE=LISINC(/2)
  516.          NLIGRF=LISDUA(/2)
  517.          SEGINI,INCTRA
  518.  
  519.          INCTRR(IRI)=INCTRA
  520.  
  521.          MELEME=IRIGEL(1,IRI)
  522.          SEGACT,MELEME
  523.          N2=NUM(/2)
  524.  
  525.          xmatri = irigel(4,iri)
  526.          segact xmatri
  527.  
  528.          DO 34 J=1,NLIGRE
  529.             DO 33 K=1,MAXINC
  530.                IF(LISINC(J).NE.IMIK(K)) GO TO 33
  531.                IF(IHAR(K).EQ.IHARIR) GO TO 32
  532.    33       CONTINUE
  533.    32       CONTINUE
  534. C       ... K est tel que LISINC(J)=IMIK(K) et IHARIR=IHAR(K),
  535. C           on le met dans INCTRA(J) (J numérote les variables) correspondant ...
  536.             INCTRA(J)=K
  537. C       ... Dans la boucle ci-dessous on met à 1 les INCPO correspondants à la
  538. C           variable K pour les noeuds des éléments du maillage ...
  539.             DO 31 I=1,N2
  540.                IJ=INUINV(NUM(NOELEP(J),I))
  541.                INCPO(K,IJ)=1
  542. *  terme diagonal
  543.         if (j.le.nligrf) diatmp(K,IJ)=diatmp(k,ij)+
  544.      >    re(j,j,i)*coerig(iri)
  545.  
  546.    31       continue
  547.    34    CONTINUE
  548.  
  549.          SEGDES,INCTRA
  550.  
  551.          NLIGRF=LISINC(/2)
  552.          NLIGRE=LISDUA(/2)
  553.          SEGINI,INCTRA
  554.          INCTRD(IRI)=INCTRA
  555.  
  556.          DO 342 J=1,NLIGRE
  557.             DO 332 K=1,MAXDUA
  558.                IF(LISDUA(J).NE.IDUA(K)) GO TO 332
  559.                IF(MHAR1.IHAR(K).EQ.IHARIR) GO TO 322
  560.   332       CONTINUE
  561.   322       CONTINUE
  562. C       ... K est tel que LISDUA(J)=IDUA(K) et IHARIR=IHAR(K),
  563. C           on le met dans INCTRA(J) (J numérote les variables) correspondant ...
  564.             INCTRA(J)=K
  565. C       ... Dans la boucle ci-dessous on met à 1 les INCPO correspondants à la
  566. C           variable K pour les noeuds des éléments du maillage ...
  567.             DO  I=1,N2
  568.                IJ=INUINV(NUM(NOELED(J),I))
  569.                MIPO1.INCPO(K,IJ)=1
  570. *  terme diagonal
  571.                if (j.le.nligrf) diatmp(K,IJ)=diatmp(k,ij)+
  572.      >    re(j,j,i)*coerig(iri)
  573.             enddo
  574.   342    CONTINUE
  575.  
  576.          SEGDES,DESCR
  577.          SEGDES,INCTRA
  578.          SEGDES,MELEME
  579.    29 CONTINUE
  580. C
  581. C  ****  INITIALISATION DE IPOS
  582. C
  583. C ... IPOS(I+1)-IPOS(I) = nombre de colonnes liées au noeud I ...
  584. C ... IPOS(I)+1 = numéro de la première colonne concernant le noeud I ...
  585.       IPOS(1)=0
  586.       IPOD(1)=0
  587. C ... NA = nombre de 1 dans INCPO => nombre de colonnes de la matrice ...
  588.       NA=0
  589.       ND=0
  590.       DO 37 I=1,NNOE
  591.       nad=na
  592.       ndd=nd
  593.       diamax=0.d0
  594.       DO 35 K=1,MAXINC
  595.       IF(INCPO(K,I).EQ.0) GO TO 35
  596.       NA=NA+1
  597.       INCPO(K,I)=NA
  598.       itrv1(na-nad)=k  
  599.       dtrv1(na-nad)= -diatmp(k,i)
  600.       diamax=max(diamax,abs(dtrv1(na-nad)))
  601.    35 CONTINUE
  602.       diaref = diamax * xszpre
  603.       do k=1,na-nad
  604.        if (abs(dtrv1(k)).lt.diaref) then
  605. **      write (6,*) ' terme diag petit ',dtrv1(k)
  606.         dtrv1(k)=dtrv1(k)+diamax
  607.        endif
  608.       enddo
  609. *  trier incpo suivant les val de diatmp
  610.       call triflo(dtrv1,dtrv2,itrv1,itrv2,na-nad)
  611.       do k=1,na-nad
  612.        incpo(itrv1(k),i)=k+nad
  613.       enddo
  614.       IPOS(I+1)=NA
  615. C ... ND = nombre de 1 dans MIPO1.INCPO => nombre de lignes de la matrice ...
  616.       diamax=0.d0
  617.          DO 352 K=1,MAXDUA
  618.             IF(MIPO1.INCPO(K,I).NE.0) THEN
  619.                ND=ND+1
  620. C          ... MIPO1.INCPO(K,I) = numéro de l'équation ...
  621.                MIPO1.INCPO(K,I)=ND
  622.          itrv1(nd-ndd)=k
  623.          dtrv1(nd-ndd)= -diatmp(k,i)
  624.          diamax=max(diamax,abs(dtrv1(nd-ndd)))
  625.             ENDIF
  626.   352    CONTINUE
  627.       diaref = diamax * xszpre
  628.       do k=1,nd-ndd
  629.        if (abs(dtrv1(k)).lt.diaref) then
  630. **      write (6,*) ' terme diag petit ',dtrv1(k)
  631.         dtrv1(k)=dtrv1(k)+diamax
  632.        endif
  633.       enddo
  634.       call triflo(dtrv1,dtrv2,itrv1,itrv2,nd-ndd)
  635.       do k=1,nd-ndd
  636.        mipo1.incpo(itrv1(k),i)=k+ndd
  637.       enddo
  638.          IPOD(I+1)=ND
  639.    37 CONTINUE
  640. *        write(*,*) 'Nb de colonnes de la matrice : ',NA,maxinc
  641. *        write(*,*) 'Nb de lignes   de la matrice : ',ND,maxdua
  642.  
  643.  
  644.       SEGDES,MIDUA,MIMIK,MHARK,MHAR1
  645.  
  646. C ... On va tester que tout est OK pour la suite ...
  647.  
  648.        IF(NA.NE.ND) THEN
  649. *          write(6,*) ' ipos'
  650. *          write(6,*) ( ipos(IU),IU=1,ipos(/1))
  651. *          write(6,*) ' ipod '
  652. *          write(6,*) ( ipod(IU),IU=1,ipod(/1))
  653.           CALL ERREUR(756)
  654.           RETURN
  655.        ENDIF
  656.  
  657.        DO 567 IINO=1,NNOE1
  658.           IF(IPOS(IINO).NE.IPOD(IINO)) THEN
  659.              WRITE(*,*) 'ERREUR dans ASNS1 !!! IPOS != IPOD !!!'
  660.              RETURN
  661.           ENDIF
  662.  567  CONTINUE
  663. C
  664. C  ****  INITIALISATION DE IMINI a été supprimée car ce segment
  665. C        ne servait à rien ...
  666. *      write(6,*) ' ipos', ( ipos(iu),iu=1,ipos(/1))
  667. *      write(6,*) ' ipod', ( ipod(iu),iu=1,ipod(/1))
  668. *      write(6,*) ' itopo', ( itopo(iu),iu=1,itopo(/1))
  669. *      write(6,*) ' itopob', ( itopob(iu),iu=1,itopob(/1))
  670. *      write(6,*) ' iitop', ( iitop(iu),iu=1,iitop(/1))
  671. *      write(6,*) ' iitopb', ( iitopb(iu),iu=1,iitopb(/1))
  672.       SEGsup DIATMP,strv
  673.       SEGDES,MRIGID
  674.       SEGDES,IPOS,IPOD
  675.       SEGDES,ITOPO,ITOPOB
  676.       SEGDES,IITOP,IITOPB
  677.       SEGDES,INUINV
  678.       SEGDES,IPT1
  679.       SEGDES,MINCPO
  680.       SEGDES,MIPO1
  681.       SEGSUP,IMIN,IMINB
  682.       SEGDES,INCTRR
  683.       INCTRY=INCTRR
  684.       SEGDES,INCTRD
  685.       INCTRZ=INCTRD
  686.       SEGINI,MMATRI
  687.       NENS=0
  688.       IGEOMA=IPT1
  689.       IIDUA=MIDUA
  690.       IINCPO=MINCPO
  691.       IDUAPO=MIPO1
  692.       IIMIK=MIMIK
  693.       IHARK=MHARK
  694.       IHARDU=MHAR1
  695.       INUINY=INUINV
  696.       ITOPOY=ITOPO
  697.       ITOPOD=ITOPOB
  698.       IITOPD=IITOPB
  699.       IITOPY=IITOP
  700.       MMATRX=MMATRI
  701. ccc      IMINIY=IMINI
  702.          iminiy=0
  703.       IPOY=IPOS
  704.       IPODD=IPOD
  705.       SEGDES,MMATRI
  706.       RETURN
  707.       END
  708.  
  709.  
  710.  
  711.  
  712.  
  713.  
  714.  
  715.  
  716.  
  717.  
  718.  
  719.  
  720.  
  721.  
  722.  
  723.  
  724.  
  725.  
  726.  
  727.  
  728.  
  729.  
  730.  
  731.  
  732.  
  733.  
  734.  
  735.  
  736.  
  737.  
  738.  
  739.  
  740.  

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