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

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