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    CB215821  20/11/25    13:34:15     10792          
  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.          SEGDES KMORS
  138.          SEGINI,KISA=AISA
  139.          SEGDES KISA
  140.       ELSE
  141. C   1) On s'occupe du second membre et de l'inconnue
  142. C      On parcourt le chpoint un peu comme dans ch2vec.eso
  143. C      Les DDL ou on applique une condition de Dirichlet
  144. C      ont un DDLCL(No DDL)=1 et 0 sinon
  145.          IF (IMPR.GT.6) THEN
  146.             WRITE(IOIMP,*)
  147.      $           'Dirichlet sur l''inconnue et le second membre'
  148.          ENDIF
  149.          JG=NTTT
  150.          SEGINI DDLPD
  151.          SEGINI DDLDP
  152.          SEGACT KS2B*MOD
  153.          IF (LINCX) SEGACT INCX*MOD
  154.          MCHPOI=KCLIM
  155.          SEGACT MINC
  156.          NBI=LISINC(/2)
  157. C In KRIPAD : SEGACT ISPGT
  158. C             SEGINI MLENTI
  159.          CALL KRIPAD(ISPGT,MLENTI)
  160.          SEGDES ISPGT
  161.          SEGACT IDMATP
  162.          SEGACT IDMATD
  163.          SEGACT MCHPOI
  164.          NSOUPO=IPCHP(/1)
  165.          DO 1 ISOUPO=1,NSOUPO
  166.             MSOUPO=IPCHP(ISOUPO)
  167.             SEGACT MSOUPO
  168.             NC=NOCOMP(/2)
  169.             MELEME=IGEOC
  170.             MPOVAL=IPOVAL
  171.             SEGACT MELEME
  172.             SEGACT MPOVAL
  173. C         NPTI=VPOCHA(/1)
  174.             DO 11 INC=1,NC
  175.                NOMINC=NOCOMP(INC)
  176.                FLINC=.FALSE.
  177. C  Repeat..until
  178.                INBI=1
  179.  111           CONTINUE
  180.                IF (NOMINC.EQ.LISINC(INBI)) THEN
  181.                   FLINC=.TRUE.
  182.                ELSEIF (INBI.LT.NBI) THEN
  183.                   INBI=INBI+1
  184.                   GOTO 111
  185.                ENDIF
  186.                IF (.NOT.FLINC) THEN
  187. *                  WRITE(IOIMP,*) 'melim.eso : composante ',NOMINC
  188. *     $                 ,'unknown'
  189.                   WRITE(IOIMP,*) 'CLIM : component name ',NOMINC
  190.      $                 ,' not in the matrix'
  191.                   GOTO 9999
  192.                ELSE
  193.                   N=VPOCHA(/1)
  194.                   DO 112 IN=1,N
  195.                      I1=LECT(NUM(1,IN))
  196.                      IF(I1.EQ.0)THEN
  197. *                        WRITE(IOIMP,*) 'melim.eso : le point ',IN
  198. *     $                       ,' du ch.'
  199. *                        WRITE(IOIMP,*) 'n''appartient pas au vec.'
  200.                   WRITE(IOIMP,*) 'CLIM : point number ',(num(1,IN)),
  201.      $                       ' unknown ',NOMINC,' not in the matrix'
  202.                         GOTO 9999
  203.                      ELSE
  204.                         IF (MPOS(I1,INBI).EQ.0) THEN
  205.                      WRITE(IOIMP,*) 'CLIM : point number ',(num(1,IN)),
  206.      $                          ' unknown ',NOMINC,' not in the matrix'
  207.                            GOTO 9999
  208.                         ELSE
  209.                            IDIAGP=IDMATP.NUAN(NPOS(I1)+MPOS(I1,INBI)-1)
  210.                            IDIAGD=IDMATD.NUAN(NPOS(I1)+MPOS(I1,INBI)-1)
  211. C*                        INBVA=NUAN(NPOS(I1)+MPOS(I1,INBI)-1)
  212. C                     INPT=NUAN(I1)
  213. C                     INBVA=NPOS(INPT)+MPOS(INPT,INBI)-1
  214. C                    KVEC.A(INBVA)=KVEC.A(INBVA)+VPOCHA(IN,INC)
  215. C Je préfère surcharger
  216. C
  217. *                        KS2B.A(INBVA)=VPOCHA(IN,INC)
  218. *                        IF (LINCX) INCX.A(INBVA)=VPOCHA(IN,INC)
  219.                            KS2B.A(IDIAGD)=VPOCHA(IN,INC)
  220.                            IF (LINCX) INCX.A(IDIAGP)=VPOCHA(IN,INC)
  221.                            DDLPD.LECT(IDIAGP)=IDIAGD
  222.                            DDLDP.LECT(IDIAGD)=IDIAGP
  223. *                        DDLCL.LOGI(INBVA)=.TRUE.
  224.                         ENDIF
  225.                      ENDIF
  226.  112              CONTINUE
  227.                ENDIF
  228.  11         CONTINUE
  229.             SEGDES MSOUPO
  230.             SEGDES MELEME
  231.             SEGDES MPOVAL
  232.  1       CONTINUE
  233.          SEGDES MCHPOI
  234.          SEGDES IDMATD
  235.          SEGDES IDMATP
  236.          SEGSUP MLENTI
  237.          SEGDES MINC
  238.          IF (IMPR.GT.7) THEN
  239.             SEGPRT,KS2B
  240.             IF (LINCX) SEGPRT,INCX
  241.          ENDIF
  242.          IF (LINCX) SEGDES,INCX
  243. C   2) On s'occupe de la matrice :
  244. C      pour les no de DDL ou on impose une CL dirichlet, il faut :
  245. C      annuler la ligne correspondante sauf la diagonale qui vaut 1
  246. C      reporter -colonne*clim sur le second membre
  247. C      On parcourt donc toute la matrice
  248.          IF (IMPR.GT.6) THEN
  249.             WRITE(IOIMP,*) 'Dirichlet sur la matrice'
  250.          ENDIF
  251.          SEGINI,KMORS=AMORS
  252.          SEGINI,KISA=AISA
  253.          DO 2 INTTT=1,NTTT
  254. *!            IF (DDLCL.LOGI(INTTT)) THEN
  255.             IDIAGP=DDLDP.LECT(INTTT)
  256.             IF (IDIAGP.NE.0) THEN
  257. C
  258. C     CL de dirichlet pour la ligne INTTT :
  259. C     On met 1 sur la diagonale et 0 ailleurs,
  260. C     le second membre ne change pas
  261. C
  262.                FLDIA=.FALSE.
  263.                DO 21 ICOL=KMORS.IA(INTTT),(KMORS.IA(INTTT+1)-1)
  264. *!                  IF (KMORS.JA(ICOL).NE.INTTT) THEN
  265.                   IF (KMORS.JA(ICOL).NE.IDIAGP) THEN
  266.                      KISA.A(ICOL)=ZERO
  267.                   ELSE
  268.                      KISA.A(ICOL)=ONE
  269.                      FLDIA=.TRUE.
  270.                   ENDIF
  271.  21            CONTINUE
  272.                IF (.NOT.FLDIA) THEN
  273. *!                  WRITE(IOIMP,*) 'melim.eso : diag.',INTTT,'inexistante'
  274.                   WRITE(IOIMP,*) 'melim.eso : diag.',IDIAGP,
  275.      $                 'inexistante'
  276.                   GOTO 9999
  277.                ENDIF
  278.             ELSE
  279. C
  280. C     Ligne INTTT sans condition de dirichlet :
  281. C     On annule les termes des colonnes liées à une CL dirichlet,
  282. C     et on reporte le produit colonne*CL au second membre
  283. C
  284.                TEMP=ZERO
  285.                DO 22 ICOL=KMORS.IA(INTTT),(KMORS.IA(INTTT+1)-1)
  286.                   JNTTT=KMORS.JA(ICOL)
  287.                   IDIAGD=DDLPD.LECT(JNTTT)
  288. *!                  IF (DDLCL.LOGI(JNTTT)) THEN
  289.                   IF (IDIAGD.NE.0) THEN
  290. *!                     IF (KS2B.A(JNTTT).NE.ZERO) THEN
  291. *!                     TEMP=TEMP+(KISA.A(ICOL)*KS2B.A(JNTTT))
  292.                      TEMP=TEMP+(KISA.A(ICOL)*KS2B.A(IDIAGD))
  293. *!                     ENDIF
  294.                      KISA.A(ICOL)=ZERO
  295.                   ENDIF
  296.  22            CONTINUE
  297.                IF (TEMP.NE.ZERO) THEN
  298.                   KS2B.A(INTTT)=KS2B.A(INTTT)-TEMP
  299.                ENDIF
  300.             ENDIF
  301.  2       CONTINUE
  302.          SEGDES KISA
  303.          SEGDES KMORS
  304.          SEGDES KS2B
  305.          SEGSUP DDLPD
  306.          SEGSUP DDLDP
  307.       ENDIF
  308. *
  309. *     Terminaison normale
  310. *
  311.       IRET=0
  312.       RETURN
  313. *
  314. * Format handling
  315. *
  316.  1002 FORMAT(10(1X,1PE11.4))
  317. *
  318. * Error handling
  319. *
  320.  9999 CONTINUE
  321.       IRET=1
  322.       WRITE(IOIMP,*) 'An error was detected in melim.eso'
  323.       RETURN
  324. *
  325. *     End of MELIM
  326. *
  327.       END
  328.  
  329.  

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