Télécharger la plaquette Cast3M


fr en

Télécharger melim.eso

Retour à la liste

Numérotation des lignes : 

melim

  1. C MELIM     SOURCE    GOUNAND   25/04/30    21:15:20     12258          
  2.       SUBROUTINE MELIM(MATRIK,KCLIM,
  3.      $     INCX,KS2B,
  4.      $     KMORS,KISA,
  5.      $     IMPR,IRET)
  6.       IMPLICIT REAL*8 (A-H,O-Z)
  7.       IMPLICIT INTEGER (I-N)
  8. C***********************************************************************
  9. C NOM         : MELIM
  10. C DESCRIPTION :
  11. C     Ce SP s'occupe de prendre en compte les  conditions
  12. C     aux limites de type Dirichlet pour un système linéaire
  13. C     du type Ax=b (avec A utilisant le stockage Morse.
  14. C     On procède en trois étapes :
  15. C     1 - on surcharge le second membre (b)
  16. C         et l'inconnue (x) avec le chpoint des conditions
  17. C         aux limites
  18. C     2 - on s'occupe de la matrice Morse (A)
  19. C         Aux degrés de libertés où il y a des
  20. C         conditions aux limites     : 1 sur la diagonale
  21. C                                      0 sur la ligne
  22. C                                      0 sur la colonne en reportant les
  23. C                                        produits dans le second membre
  24. C                                        (sinon perte de symétrie
  25. C                                         éventuelle de A)
  26. C     3 - On nettoie la matrice morse A des 0 qu'on vient de lui
  27. C         rajouter (sous-programme clmors).
  28. C
  29. C     2 et 3 peuvent sans doute etre regroupées au prix d'une
  30. C            perte de lisibilité
  31. C
  32. C
  33. C LANGAGE     : ESOPE
  34. C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF)
  35. C               mél : gounand@semt2.smts.cea.fr
  36. C***********************************************************************
  37. C APPELES  : KRIPAD, CLMORS
  38. C***********************************************************************
  39. C ENTREES            : MATRIK, KCLIM, IMPR
  40. C ENTREES/SORTIES    : KS2B, INCX
  41. C SORTIES            : KMORS,KISA,IRET
  42. C CODE RETOUR (IRET) : 0 si ok
  43. C                     <0 si problème
  44. C     MATRIK : pointeur sur segment MATRIK de l'include SMMATRIK
  45. C              on pioche dedans les informations nécessaires
  46. C              (différents pointeurs, nb. de ddl...)
  47. C     KCLIM  : pointeur sur segment MCHPOI de l'include SMCHPOI
  48. C              chpoint contenant les conditions aux limites de
  49. C              Dirichlet (valeur imposée).
  50. C     IMPR   : niveau d'impression
  51. C     KS2B   : pointeur sur segment IZA de l'include SMMATRIK
  52. C              vecteur second membre
  53. C              (ie b dans la résolution de Ax=b).
  54. C     INCX   : pointeur sur segment IZA de l'include SMMATRIK
  55. C              vecteur des inconnues
  56. C              (ie x dans la résolution de Ax=b).
  57. C     KMORS  : pointeurs sur segment PMORS et IZA de l'include
  58. C     KISA     SMMATRIK : profil et valeurs de la matrice morse
  59. C              (ie A dans la résolution de Ax=b).
  60. C***********************************************************************
  61. C VERSION    : v1, 20/12/99
  62. C HISTORIQUE : v1, 01/04/98, création
  63. C HISTORIQUE :     20/12/99,
  64. C   Nouvelle signification de NUAN (renumérotation des ddl et non pas
  65. C   des points)
  66. C HISTORIQUE : 05/01/00 : non-duplication de la matrice assemblée
  67. C HISTORIQUE : 09/04/04 : idmatp idmatd,
  68. C                         la cl dirichlet n'est plus forcément sur la
  69. C                         diagonale
  70. C HISTORIQUE :
  71. C***********************************************************************
  72. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  73. C en cas de modification de ce sous-programme afin de faciliter
  74. C la maintenance !
  75. C***********************************************************************
  76.  
  77. -INC PPARAM
  78. -INC CCOPTIO
  79. -INC SMELEME
  80.       POINTEUR ISPGT.MELEME
  81. -INC SMLENTI
  82.       POINTEUR DDLPD.MLENTI
  83.       POINTEUR DDLDP.MLENTI
  84. *-INC SMLLOGI
  85. *      SEGMENT MLLOGI
  86. *      LOGICAL LOGI(JG)
  87. *      ENDSEGMENT
  88. *      INTEGER JG
  89. *      POINTEUR DDLCL.MLLOGI
  90. -INC SMCHPOI
  91.       POINTEUR KCLIM.MCHPOI
  93.       POINTEUR KMORS.PMORS
  94.       POINTEUR KISA.IZA
  95.       POINTEUR AMORS.PMORS
  96.       POINTEUR AISA.IZA
  97.       POINTEUR KS2B.IZA
  98.       POINTEUR INCX.IZA
  99.       POINTEUR IDMATP.IDMAT
  100.       POINTEUR IDMATD.IDMAT
  101. *
  102.       INTEGER IMPR,IRET
  103. *
  104. *     .. Variables locales
  105. *     .. Parameters
  106.       REAL*8 ZERO,ONE
  107.       PARAMETER (ZERO=0.0D0, ONE=1.0D0)
  108. *     ..
  109.       CHARACTER*(LOCOMP) NOMINC
  110.       LOGICAL FLINC
  111.       LOGICAL FLDIA
  112. *
  113.       INTEGER I1,IN,ISOUPO,INC,INBI,INBVA
  114.       INTEGER    N ,NSOUPO,NC ,NBI
  115.       INTEGER ICOL,INTTT,JNTTT
  116.       INTEGER      NTTT
  117.       REAL*8 TEMP
  118.       LOGICAL LINCX
  119. C***
  120.       LINCX=(INCX.NE.0)
  121. C     On récupère les segments utiles
  122.       IF (IMPR.GT.5) THEN
  123.          WRITE(IOIMP,*) 'melim.eso : CL de Dirichlet'
  124.       ENDIF
  125.       SEGACT MATRIK
  126.       MINC =KMINC
  127.       ISPGT=KISPGT
  128.       NTTT =KNTTT
  129.       IDMATP=KIDMAT(1)
  130.       IDMATD=KIDMAT(2)
  131.       AMORS=KIDMAT(4)
  132.       AISA =KIDMAT(5)
  133.       SEGDES MATRIK
  134. C     S'il n'y a pas de CL de Dirichlet, on n'a rien à faire
  135.       IF(KCLIM.EQ.0) THEN
  136.          SEGINI,KMORS=AMORS
  137.          SEGINI,KISA=AISA
  138.       ELSE
  139. C   1) On s'occupe du second membre et de l'inconnue
  140. C      On parcourt le chpoint un peu comme dans ch2vec.eso
  141. C      Les DDL ou on applique une condition de Dirichlet
  142. C      ont un DDLCL(No DDL)=1 et 0 sinon
  143.          IF (IMPR.GT.6) THEN
  144.             WRITE(IOIMP,*)
  145.      $           'Dirichlet sur l''inconnue et le second membre'
  146.          ENDIF
  147.          JG=NTTT
  148.          SEGINI DDLPD
  149.          SEGINI DDLDP
  150.          MCHPOI=KCLIM
  151.          SEGACT MINC
  152.          NBI=LISINC(/2)
  153. C In KRIPAD : SEGACT ISPGT
  154. C             SEGINI MLENTI
  155.          CALL KRIPAD(ISPGT,MLENTI)
  156.          SEGDES ISPGT
  157.          SEGACT IDMATP
  158.          SEGACT IDMATD
  159.          SEGACT MCHPOI
  160.          NSOUPO=IPCHP(/1)
  161.          DO 1 ISOUPO=1,NSOUPO
  162.             MSOUPO=IPCHP(ISOUPO)
  163.             SEGACT MSOUPO
  164.             NC=NOCOMP(/2)
  165.             MELEME=IGEOC
  166.             MPOVAL=IPOVAL
  167.             SEGACT MELEME
  168.             SEGACT MPOVAL
  169. C         NPTI=VPOCHA(/1)
  170.             DO 11 INC=1,NC
  171.                NOMINC=NOCOMP(INC)
  172.                FLINC=.FALSE.
  173. C  Repeat..until
  174.                INBI=1
  175.  111           CONTINUE
  176.                IF (NOMINC.EQ.LISINC(INBI)) THEN
  177.                   FLINC=.TRUE.
  178.                ELSEIF (INBI.LT.NBI) THEN
  179.                   INBI=INBI+1
  180.                   GOTO 111
  181.                ENDIF
  182.                IF (.NOT.FLINC) THEN
  183. *                  WRITE(IOIMP,*) 'melim.eso : composante ',NOMINC
  184. *     $                 ,'unknown'
  185.                   WRITE(IOIMP,*) 'CLIM : component name ',NOMINC
  186.      $                 ,' not in the matrix'
  187.                   GOTO 9999
  188.                ELSE
  189.                   N=VPOCHA(/1)
  190.                   DO 112 IN=1,N
  191.                      I1=LECT(NUM(1,IN))
  192.                      IF(I1.EQ.0)THEN
  193. *                        WRITE(IOIMP,*) 'melim.eso : le point ',IN
  194. *     $                       ,' du ch.'
  195. *                        WRITE(IOIMP,*) 'n''appartient pas au vec.'
  196.                   WRITE(IOIMP,*) 'CLIM : point number ',(num(1,IN)),
  197.      $                       ' unknown ',NOMINC,' not in the matrix'
  198.                         GOTO 9999
  199.                      ELSE
  200.                         IF (MPOS(I1,INBI).EQ.0) THEN
  201.                      WRITE(IOIMP,*) 'CLIM : point number ',(num(1,IN)),
  202.      $                          ' unknown ',NOMINC,' not in the matrix'
  203.                            GOTO 9999
  204.                         ELSE
  205.                            IDIAGP=IDMATP.NUAN(NPOS(I1)+MPOS(I1,INBI)-1)
  206.                            IDIAGD=IDMATD.NUAN(NPOS(I1)+MPOS(I1,INBI)-1)
  207. C*                        INBVA=NUAN(NPOS(I1)+MPOS(I1,INBI)-1)
  208. C                     INPT=NUAN(I1)
  209. C                     INBVA=NPOS(INPT)+MPOS(INPT,INBI)-1
  210. C                    KVEC.A(INBVA)=KVEC.A(INBVA)+VPOCHA(IN,INC)
  211. C Je préfère surcharger
  212. C
  213. *                        KS2B.A(INBVA)=VPOCHA(IN,INC)
  214. *                        IF (LINCX) INCX.A(INBVA)=VPOCHA(IN,INC)
  215.                            KS2B.A(IDIAGD)=VPOCHA(IN,INC)
  216.                            IF (LINCX) INCX.A(IDIAGP)=VPOCHA(IN,INC)
  217.                            DDLPD.LECT(IDIAGP)=IDIAGD
  218.                            DDLDP.LECT(IDIAGD)=IDIAGP
  219. *                        DDLCL.LOGI(INBVA)=.TRUE.
  220.                         ENDIF
  221.                      ENDIF
  222.  112              CONTINUE
  223.                ENDIF
  224.  11         CONTINUE
  225.             SEGDES MSOUPO
  226.             SEGDES MELEME
  227.             SEGDES MPOVAL
  228.  1       CONTINUE
  229.          SEGDES MCHPOI
  230.          SEGDES IDMATD
  231.          SEGDES IDMATP
  232.          SEGSUP MLENTI
  233.          SEGDES MINC
  234.          IF (IMPR.GT.7) THEN
  235.             SEGPRT,KS2B
  236.             IF (LINCX) SEGPRT,INCX
  237.          ENDIF
  238. C   2) On s'occupe de la matrice :
  239. C      pour les no de DDL ou on impose une CL dirichlet, il faut :
  240. C      annuler la ligne correspondante sauf la diagonale qui vaut 1
  241. C      reporter -colonne*clim sur le second membre
  242. C      On parcourt donc toute la matrice
  243.          IF (IMPR.GT.6) THEN
  244.             WRITE(IOIMP,*) 'Dirichlet sur la matrice'
  245.          ENDIF
  246.          SEGINI,KMORS=AMORS
  247.          SEGINI,KISA=AISA
  248.          DO 2 INTTT=1,NTTT
  249. *!            IF (DDLCL.LOGI(INTTT)) THEN
  250.             IDIAGP=DDLDP.LECT(INTTT)
  251.             IF (IDIAGP.NE.0) THEN
  252. C
  253. C     CL de dirichlet pour la ligne INTTT :
  254. C     On met 1 sur la diagonale et 0 ailleurs,
  255. C     le second membre ne change pas
  256. C
  257.                FLDIA=.FALSE.
  258.                DO 21 ICOL=KMORS.IA(INTTT),(KMORS.IA(INTTT+1)-1)
  259. *!                  IF (KMORS.JA(ICOL).NE.INTTT) THEN
  260.                   IF (KMORS.JA(ICOL).NE.IDIAGP) THEN
  261.                      KISA.A(ICOL)=ZERO
  262.                   ELSE
  263.                      KISA.A(ICOL)=ONE
  264.                      FLDIA=.TRUE.
  265.                   ENDIF
  266.  21            CONTINUE
  267.                IF (.NOT.FLDIA) THEN
  268. *!                  WRITE(IOIMP,*) 'melim.eso : diag.',INTTT,'inexistante'
  269.                   WRITE(IOIMP,*) 'melim.eso : diag.',IDIAGP,
  270.      $                 'inexistante'
  271.                   GOTO 9999
  272.                ENDIF
  273.             ELSE
  274. C
  275. C     Ligne INTTT sans condition de dirichlet :
  276. C     On annule les termes des colonnes liées à une CL dirichlet,
  277. C     et on reporte le produit colonne*CL au second membre
  278. C
  279.                TEMP=ZERO
  280.                DO 22 ICOL=KMORS.IA(INTTT),(KMORS.IA(INTTT+1)-1)
  281.                   JNTTT=KMORS.JA(ICOL)
  282.                   IDIAGD=DDLPD.LECT(JNTTT)
  283. *!                  IF (DDLCL.LOGI(JNTTT)) THEN
  284.                   IF (IDIAGD.NE.0) THEN
  285. *!                     IF (KS2B.A(JNTTT).NE.ZERO) THEN
  286. *!                     TEMP=TEMP+(KISA.A(ICOL)*KS2B.A(JNTTT))
  287.                      TEMP=TEMP+(KISA.A(ICOL)*KS2B.A(IDIAGD))
  288. *!                     ENDIF
  289.                      KISA.A(ICOL)=ZERO
  290.                   ENDIF
  291.  22            CONTINUE
  292.                IF (TEMP.NE.ZERO) THEN
  293.                   KS2B.A(INTTT)=KS2B.A(INTTT)-TEMP
  294.                ENDIF
  295.             ENDIF
  296.  2       CONTINUE
  297.          SEGSUP DDLPD
  298.          SEGSUP DDLDP
  299.       ENDIF
  300. *
  301. *     Terminaison normale
  302. *
  303.       IRET=0
  304.       RETURN
  305. *
  306. * Format handling
  307. *
  308.  1002 FORMAT(10(1X,1PE11.4))
  309. *
  310. * Error handling
  311. *
  312.  9999 CONTINUE
  313.       IRET=1
  314.       WRITE(IOIMP,*) 'An error was detected in melim.eso'
  315.       RETURN
  316. *
  317. *     End of MELIM
  318. *
  319.       END
  320.  
  321.  

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