Télécharger la plaquette Cast3M


fr en

Télécharger profch.eso

Retour à la liste

Numérotation des lignes : 

profch

  1. C PROFCH    SOURCE    CB215821  26/08/24    21:17:58     12622          
  2.       SUBROUTINE PROFCH(MATRAK,IMPR,KPIMP,NIMP)
  3. C************************************************************************
  4. C
  5. C     OBJET :
  6. C        CE SOUS PROGRAMME PREPARE L'INVERSION DE LA MATRICE DE PRESSION
  7. C
  8. C        KTYPI= 1 OU 5, METHODE DIRECTE
  9. C
  10. C                       1/ IOPTI = 0 PAS DE RENUMEROTATION
  11. C                          IOPTI = 1 RENUMEROTATION DE LA MATRICE
  12. C                                    DE PRESSION
  13. C                       2/ CALCUL DU PROFIL ET DE LA TAILLE
  14. C                          DE LA DEMI-MATRICE
  15. C
  16. C        MATRAK : objet de type MATRAK
  17. C                 Si KIZCL=0 On cree IZL
  18. C                 Sinon rien
  19. C        IMPR 1 ou 2 niveau d'impressions
  20. C        KPIMP = 0 pas de pression imposee
  21. C        KPIMP = 1 pression imposée en un point (le dernier)
  22. C        KPIMP = 2 pression moyenne imposée
  23. C        NIMP      numero de la pression imposée
  24. C
  25. C************************************************************************
  26.       IMPLICIT INTEGER(I-N)
  27.       IMPLICIT REAL*8 (A-H,O-Z)
  28.  
  29. C***
  30. C-INC SMMATRAKANC
  31. C*************************************************************************
  32. C
  33. C    REPERAGE ET STOKAGE DES MATRICES ELEMENTAIRES  puis assemblees
  34. C
  35.  
  36. * LGEOC SPG de la pression et/ou des multiplicateurs de Lagrange
  37. * (points CENTRE ) pour chaque operateur de contrainte
  38. * KGEOC SPG pour la totalite des points CENTRE.
  39. * KGEOS SPG pour la totalite des points SOMMET (Diagonale vitesse)
  40. * KLEMC Connectivites de l'ensemble des contraintes
  41. * LIZAFM(NBSOUS) contient les pointeurs IZAFM des sous-zones
  42.  
  43.       SEGMENT MATRAK
  44.       INTEGER LGEOC(NBOP),IDEBS(NBOP),IFINS(NBOP)
  45.       INTEGER LIZAFM(NBSOUS)
  46.       INTEGER IKAM0 (NBSOUS)
  47.       INTEGER IMEM  (NBELC)
  48.       INTEGER KLEMC,KGEOS,KGEOC,KDIAG,KCAC,KIZCL,KIZGC
  49.       ENDSEGMENT
  50.  
  51.       SEGMENT IZAFM
  52.       REAL*8  AM(NNELP,NP,IESP),RPGI(NELAX)
  53.       ENDSEGMENT
  54.  
  55.       POINTEUR IPMJ.IZAFM,IPMK.IZAFM
  56.  
  57. C*******************************************************************
  59. -INC SMCOORD
  60. -INC SMLENTI
  61.       POINTEUR IZIPAD.MLENTI
  62. -INC SMELEME
  63.       POINTEUR IZK.MELEME
  64. -INC SMTABLE
  65.       POINTEUR MTABP.MTABLE
  66.       PARAMETER (IBLSIZ=50000)
  67.  
  68.       SEGMENT IZTRAV
  69.          INTEGER ITAB(NJA)
  70.       ENDSEGMENT
  71. C
  72. C IVK = Tableau de travail
  73.       PARAMETER (MAXIVK=100)
  74.       INTEGER IVK(MAXIVK)
  75. C***
  76.       CHARACTER*8 TYPE
  77.       CHARACTER*4 LMOT(1)
  78.       DATA LMOT/'IMPR'/
  79.       DATA IOPTI/1/
  80. C***
  81.  
  82.       SEGACT MATRAK*MOD
  83.       IF(KIZCL.NE.0)THEN
  84.       SEGDES MATRAK
  85.       RETURN
  86.       ENDIF
  87.  
  88.       KGEO=KGEOC
  89.       CALL KNBEL(KGEO,NBELC)
  90.       MELEME=KLEMC
  91.  
  92. C
  93. C--- CALNOE renvoie dans IZTKB un tableau contenant pour chaque noeud
  94. C    du domaine la liste des {l{ments poss{dant ce noeud. Les noeuds
  95. C    sont rep{r{s dans la num{rotation locale(domaniale), les {l{ments
  96. C    sont pris dans l'ordre fourni par la num{rotation naturelle.
  97. C    Il est possible, partant de ce tableau, de construire pour chaque
  98. C    {l{ment la liste des {l{ments voisins. Ceci est tr}s utile pour
  99. C    le calcul de la matrice de pression ou Akl.ne.0 si k et l sont
  100. C    deux {l{ments voisins.
  101. C
  102. C    Cette technique n'est pas (encore) utilis{e si KTYPI = 1 ou 5.
  103. C    Les maillages trait{s dans ce cas {tant "petits", ce n'est pas
  104. C    trop grave. Par contre lorsque KTYPI=2,3,4, c'est @ dire pour les
  105. C    tr}s gros maillages, son emploi est obligatoire.
  106. C
  107. C    Cet appel ne doit etre fait que si on stocke la matrice de
  108. C    pression (KTPI.LE.5).
  109. C
  110.       CALL KALNO0(MELEME,IZK,MLENTI,IZIPAD,IRET)
  111.       IF(IRET.EQ.0)GO TO 90
  112.  
  113. C---  Activation de tous les MELEMEs du domaine
  114.  
  115.       SEGACT MELEME,IZK,MLENTI
  116.       NBSOUS=LISOUS(/1)
  117.       IF(NBSOUS.EQ.0)NBSOUS=1
  118.       DO 400 NS=1,NBSOUS
  119.       IF(NBSOUS.EQ.1)IPT1=MELEME
  120.       IF(NBSOUS.NE.1)IPT1=LISOUS(NS)
  121.       SEGACT IPT1
  122.  400  CONTINUE
  123.  
  124. C
  125. C************************************************************************
  126. C     C H O L E S K I
  127. C************************************************************************
  128. C
  129.       NJA=NBELC
  130.       NJAN=NBELC
  131.       NJAB=NBELC
  132.  
  133.       SEGINI IZL,IZTRAV
  134.       KIZCL=IZL
  135.  
  136.       LONG=0
  137.  
  138.       CALL INITI(NUAN,NBELC,0)
  139.       CALL INITI(NUNA,NBELC,0)
  140.       NUAN(1)=1
  141.       NUNA(1)=1
  142.  
  143. C********* ON STOKE LES POINTEURS DES MELEME
  144.  
  145.       KI=1
  146.       DO 40 NS=1,NBSOUS
  147.       IF(NBSOUS.EQ.1)IPT1=MELEME
  148.       IF(NBSOUS.NE.1)IPT1=LISOUS(NS)
  149.  
  150.       NEL=IPT1.NUM(/2)
  151.       IPOINT=IPT1
  152.       CALL INITI(IMEL(KI),NEL,IPOINT)
  153.  
  154.       DO 41 I=1,NEL
  155.       IMJ(KI+I-1)=I
  156.  41   CONTINUE
  157.  
  158.       KI=KI+NEL
  159.  
  160.  40   CONTINUE
  161.  
  162. C*********
  163.  
  164.       IF(IOPTI.EQ.1)GO TO 33
  165.  
  166. C********* ON NE RENUMEROTE PAS
  167.  
  168.       IF(IMPR.GE.2)WRITE(6,*)' PAS DE RENUMEROTATION DE LA PRESSION '
  169.       DO 32 I=1,NBELC
  170.       NUAN(I)=I
  171.       NUNA(I)=I
  172.  32   CONTINUE
  173.       GO TO 34
  174.  
  175. C********* OPTIMISATION DE LA NUMEROTAION DE LA MATRICE DE PRESSION
  176.  33   CONTINUE
  177.       IF(IMPR.GE.2)WRITE(6,*)' ON RENUMEROTE LA MATRICE DE PRESSION '
  178. C*********
  179.  
  180. C     write(6,*)' SUB PROFCH : IPADL :'
  181. C     nbp=ipadl(/1)
  182. C     write(6,1001)(IPADL(II),II=1,NBP)
  183. C     im2=lect(/1)
  184. C     im1=izk.num(/1)
  185. C     write(6,*)' IM1=',IM1 ,' Taille IZK.NUM'
  186. C     do 1873 i=1,im2
  187. C     nbb=lect(i)
  188. C     write(6,*)' I=->',i,' nb=',lect(i)
  189. C     write(6,1001)(izk.num(ii,i),ii=1,nbb)
  190. C1873 continue
  191.  
  192.       KC=1
  193.  
  194.       DO 30 KK=1,NBELC
  195.  
  196.       K=NUNA(KK)
  197.  
  198.       IF(K.NE.0)GO TO 310
  199. C   ON A PLUS D'ELEMENT NOUVELLEMENT NUMEROTE
  200.       WRITE(6,*)' EST CE LA DISCONTINUITE ',K,KC,KK
  201.       KC=KC+1
  202.       DO 311 KKA=1,NBELC
  203.       IF(NUAN(KKA).EQ.0)GO TO 312
  204. C SI NUAN(KKA)=0 C'EST QUE L'ELEMENT KKA N'A PAS ETE NUMEROTE
  205.  311  CONTINUE
  206.       WRITE(6,*)' IL YA UN PROBLEME DE NUMEROTATION'
  207.       CALL ARRET(0)
  208.  312  CONTINUE
  209.       NUAN(KKA)=KC
  210.       NUNA(KC)=KKA
  211.       K=NUNA(KC)
  212.  310  CONTINUE
  213.  
  214.       IPT1=IMEL(K)
  215.       NP=IPT1.NUM(/1)
  216.       DO 1931 I=1,NP
  217.       NU=IPT1.NUM(I,IMJ(K))
  218. C     NU=IPT1.NUM(I,K)
  219.       NU=IZIPAD.LECT(NU)
  220.       NELV=LECT(NU)
  221.       DO 31 NN=1,NELV
  222.       KN=IZK.NUM(NN,NU)
  223.  
  224.       IF(NUAN(KN).NE.0)GO TO 31
  225.       KC=KC+1
  226.       NUAN(KN)=KC
  227.       NUNA(KC)=KN
  228.  31   CONTINUE
  229.  1931 CONTINUE
  230.  
  231.  30   CONTINUE
  232.  
  233.       DO 35 I=1,NBELC
  234.       ITAB(I)=IMJ(NUNA(I))
  235.  35   CONTINUE
  236.       CALL RSETI(IMJ,ITAB,NBELC)
  237.  
  238.       DO 36 I=1,NBELC
  239.       ITAB(I)=IMEL(NUNA(I))
  240.  36   CONTINUE
  241.       CALL RSETI(IMEL,ITAB,NBELC)
  242.  
  243.       SEGSUP IZTRAV
  244. C             **********************************************
  245.  34   CONTINUE
  246.  
  247.       IF(IMPR.GE.2)THEN
  248.       WRITE(6,1929)
  249.  1929 FORMAT(/10X,' SUB PROFIM : RENUMEROTAION ZONE PRESSION ')
  250.       WRITE(6,*)' CORRESPONDANCE NUMER. ANC. --> NUMER. NOUVELLE '
  251.       WRITE(6,1930)(N,NUAN(N),N=1,NBELC)
  252.       WRITE(6,*)' CORRESPONDANCE NUMER. NOUVELLE --> NUMER. ANC. '
  253.       WRITE(6,1930)(N,NUNA(N),N=1,NBELC)
  254.  1930 FORMAT(7(1X,'**',I5,3X,I5,1X,'**'))
  255.       ENDIF
  256. C
  257.       DO 10 KK=1,NBELC
  258.       K=IMJ(KK)
  259.       IPT1=IMEL(KK)
  260.       NP=IPT1.NUM(/1)
  261.       K0=KK
  262.       DO 1932 I=1,NP
  263.       NU=IPT1.NUM(I,K)
  264.       NU=IZIPAD.LECT(NU)
  265.       NELV=LECT(NU)
  266.       DO 11 NN=1,NELV
  267.       KN=IZK.NUM(NN,NU)
  268.       KN=NUAN(KN)
  269.       IF(KN.LT.K0)K0=KN
  270.  11   CONTINUE
  271.  1932 CONTINUE
  272.       IF(KK.NE.NIMP.OR.KPIMP.NE.2)THEN
  273.       KZA(KK)=KK-K0+1
  274.       ELSE
  275.       KZA(KK)=NBELC-1
  276.       ENDIF
  277.       LONG=LONG+KZA(KK)
  278.  10   CONTINUE
  279.  
  280. C
  281. C La matrice est stockée par blocs de plusieurs lignes. En effet, il faut
  282. C l'équivalent de 900 opérations pour activer/désactiver un segment. Si
  283. C une ligne de la matrice comprend 50 éléments, on fait 100 opérations
  284. C utiles et 900 inutiles lors de la résolution.
  285. C La taille des blocs est fixée par le paramètre IBLSIZ. La diagonale est
  286. C stockée dans un segment à part.
  287. C
  288. C Il nous faut constituer un certain nombre de tableaux. Le premier segment
  289. C décrit la matrice. Il contient un tableau donnant les numéros des premières
  290. C lignes de chaque bloc et un tableau de pointeurs sur des segments
  291. C descripteurs de bloc.
  292. C
  293. C Il faut d'abord déterminer le nombre de blocs.
  294. C
  295.  
  296.       NBLK=1
  297.       IF(LONG-NBELC.GT.IBLSIZ) THEN
  298.       ILBL=0
  299.       DO 12 KK=2,NBELC
  300.          LA=KZA(KK)-1
  301.          IF(LA.GT.IBLSIZ) THEN
  302.             NBLK=NBLK+2
  303.             ILBL=0
  304.             GOTO 12
  305.          ENDIF
  306.          ILBL=ILBL+LA
  307.          IF(ILBL.GT.IBLSIZ) THEN
  308.             NBLK=NBLK+1
  309.             ILBL=LA
  310.             GOTO 12
  311.          ENDIF
  312.  12   CONTINUE
  313.       ENDIF
  314. C
  315. C On peut allouer le segment descripteur de la matrice. On va pouvoir
  316. C remplir des numéros de ligne des débuts de bloc.
  317. C
  318.       SEGINI IDMAT
  319.       IF(NBLK.EQ.1) THEN
  320.          NLDBLK(1)=2
  321.          NLDBLK(2)=NBELC+1
  322.       ELSE
  323.          IBLK=1
  324.          ILBL=0
  325.          NLDBLK(1)=2
  326.          NLDBLK(NBLK+1)=NBELC+1
  327.          DO 13 KK=2,NBELC
  328.             LA=KZA(KK)-1
  329.             IF(LA.GT.IBLSIZ) THEN
  330.               NLDBLK(IBLK+1)=KK
  331.               NLDBLK(IBLK+2)=KK+1
  332.               IBLK=IBLK+2
  333.               ILBL=0
  334.               GOTO 13
  335.             ENDIF
  336.             ILBL=ILBL+LA
  337.             IF(ILBL.GT.IBLSIZ) THEN
  338.               NLDBLK(IBLK+1)=KK
  339.               ILBL=LA
  340.               IBLK=IBLK+1
  341.               GOTO 13
  342.             ENDIF
  343.  13      CONTINUE
  344.       ENDIF
  345. C
  346. C On peut maintenant remplir le segment descripteur de chaque bloc.
  347. C
  348.       DO 15 IBLK=1,NBLK
  349.       NLBLK=NLDBLK(IBLK+1)-NLDBLK(IBLK)
  350.       SEGINI IDBLK
  351.       IDESCR(IBLK)=IDBLK
  352.       ILON=0
  353.       IDEBLK(1)=1
  354.       DO 16 KK=NLDBLK(IBLK),NLDBLK(IBLK+1)-1
  355.          LA=KZA(KK)-1
  356.          ILON=ILON+LA
  357.          IDEBLK(KK-NLDBLK(IBLK)+2)=IDEBLK(KK-NLDBLK(IBLK)+1)+LA
  358.  16   CONTINUE
  359.       SEGDES IDBLK
  360.  15   CONTINUE
  361.       SEGDES IDMAT
  362. C
  363. C On stocke dans KZA1 la référence au segment IDMAT à partir duquel tout se
  364. C retrouve.
  365. C
  366. C     write(6,*)' sub PROFCH on met IDMAT -> KZA1 ',IDMAT
  367.       KZA1=IDMAT
  368.  14   CONTINUE
  369. C
  370.       IF(IMPR.GE.2)WRITE(6,*)(KZA(M),M=1,NBELC),LONG
  371. C
  372.       IF(IMPR.GE.1)WRITE(6,200)LONG
  373.       SEGDES IZL
  374.       SEGSUP IZK,IZIPAD,MLENTI
  375.  
  376. C---  D{sactivation de tous les MELEMEs
  377.  
  378.       DO 51 NS=1,NBSOUS
  379.       IF(NBSOUS.EQ.1)IPT1=MELEME
  380.       IF(NBSOUS.NE.1)IPT1=LISOUS(NS)
  381.       SEGDES IPT1
  382.  51   CONTINUE
  383.       SEGDES MELEME
  384.  
  385.       RETURN
  386.  90   CONTINUE
  387.       WRITE(6,*)' Retour anormal de PROFCH'
  388.       RETURN
  389.  
  390. C************************************************************************
  391.  100  FORMAT(40X,1H*,9X,'.... TAILLE DE LA MATRICE DE PRESSION ',I7,
  392.      *' MOTS',13X,1H*/40X,1H*,72X,1H*)
  393.  200  FORMAT(10X,1H*,9X,'.... TAILLE DE LA MATRICE DE PRESSION ',I7,
  394.      *' MOTS',13X,1H*)
  395.  1001 FORMAT(20(1X,I5))
  396.       END
  397.  
  398.  
  399.  
  400.  
  401.  
  402.  
  403.  

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