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

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