Télécharger la plaquette Cast3M


fr en

Télécharger ksprus.eso

Retour à la liste

Numérotation des lignes : 

ksprus

  1. C KSPRUS    SOURCE    CB215821  20/11/25    13:33:12     10792          
  2.       SUBROUTINE KSPRUS
  3.      &(MELEME,IPM1,IPM2,IPM3,IAXI,IKAS,INEFMD,KPRE,IZTGG1,IK1,K0,
  4.      & ANUKK,IDCEN,IKOMP,NPTU,IPADU,IZTU1,TETAN,NP,NEL,CMD)
  5.  
  6.       IMPLICIT INTEGER(I-N)
  7.       IMPLICIT REAL*8 (A-H,O-Z)
  8. C*****************************************************************************
  9. C
  10. C     Ce SP calcule les matrices elementaires de divergence alias C
  11. C
  12. C  IKAS=1 KMCT calcul de Ct  (Div U)
  13. C  IKAS=2 KMAC calcul de C uniquement (Grad p)
  14. C  IKAS=3 KCCT calcul de C assemblage pour C et Ct
  15. C
  16. C*****************************************************************************
  17.       CHARACTER*8 NOM0
  18.  
  19.  
  20. -INC PPARAM
  21. -INC CCOPTIO
  22. -INC CCGEOME
  23. -INC SMCOORD
  24. -INC SIZFFB
  25.       POINTEUR IZF1.IZFFM
  27. -INC SMLENTI
  28.       POINTEUR IPADU.MLENTI
  29. -INC SMELEME
  30. -INC SMCHPOI
  31.       POINTEUR IZTGG1.MPOVAL,IZTU1.MPOVAL,TETAN.MPOVAL
  32. -INC CCREEL
  33.       PARAMETER (LRV=64)
  34.       SEGMENT PETROV
  35.       REAL*8 WT(LRV,NP,NPG,IDIM),WS(LRV,NP,NPG,IDIM),HK(LRV,IDIM,NP,NPG)
  36.       REAL*8 UIL(LRV,IDIM,NPG),DUIL(LRV,IDIM,NPG)
  37.       REAL*8 PGSK(LRV,NPG),RPGK(LRV,NPG),AIRE(LRV),ANUK(LRV),COEFK(LRV)
  38.       REAL*8 AJK(LRV,IDIM,IDIM,NPG)
  39.       ENDSEGMENT
  40.       DIMENSION KIPM(3),XYZ1(24)
  41. C
  42. C OPERATEUR PRESSION
  43. C
  44.       DEUPI=1.D0
  45.       IF(IAXI.NE.0)DEUPI=2.D0*XPI
  46.       IF(IKAS.EQ.3)DEUPI=-DEUPI
  47.  
  48.       NBEL = NEL
  49. C     write(6,*)' DEUPI=',deupi
  50.       IF(IDIM.EQ.2)IPM3=IPM1
  51.       KIPM(1)=IPM1
  52.       KIPM(2)=IPM2
  53.       KIPM(3)=IPM3
  54.       SEGACT MELEME
  55.       SEGACT IPM1*MOD,IPM2*MOD,IPM3*MOD
  56.       SEGACT IPADU,IZTU1,TETAN,IZTGG1
  57.  
  58.       IF(KPRE.NE.2)THEN
  59. C
  60. C      ON TRAITE LE CAS DIV U -> T sommet
  61. C
  62.        IF(INEFMD.EQ.3)THEN
  63.        IF(KPRE.EQ.3)NOM0=NOMS(ITYPEL)//'PRP0'
  64.        IF(KPRE.EQ.4)NOM0=NOMS(ITYPEL)//'PRP1'
  65.        IF(KPRE.EQ.5)NOM0=NOMS(ITYPEL)//'PFP1'
  66.        ELSEIF(INEFMD.EQ.2)THEN
  67.        IF(KPRE.EQ.3)NOM0=NOMS(ITYPEL)//'MCP0'
  68.        IF(KPRE.EQ.4)NOM0=NOMS(ITYPEL)//'MCP1'
  69.        IF(KPRE.EQ.5)NOM0=NOMS(ITYPEL)//'MCF1'
  70.        ELSEIF(INEFMD.EQ.4)THEN
  71.        NOM0=NOMS(ITYPEL)//'    '
  72.        ENDIF
  73. c       write(*,*)' Kpruss : NOM0=',nom0,' ikas=',ikas,' nel=',nel
  74. c       write(*,*)' IK1',IK1,' COEF',IZTGG1.VPOCHA(1,1),'ANUK',ANUKK
  75.        CALL KALPBG(NOM0,'FONFORM ',IZFFM)
  76.  
  77.       SEGACT IZFFM*MOD
  78.       IZHR=KZHR(1)
  79.       SEGACT IZHR*MOD
  80.       NES=GR(/1)
  81.       NPG=GR(/3)
  82.       IZF1=KTP(1)
  83.       SEGACT IZF1*MOD
  84. * ON modifie pour avoir les bonnes fonctions de forme. GBM 18/10/99
  85.       MP1=FN(/1)
  86. *
  87. *     On fait des packets
  88. *
  89.       NNN=MOD(NBEL,LRV)
  90.       IF(NNN.EQ.0) THEN
  91.          NPACK=NBEL/LRV
  92.       ELSE
  93.          NPACK=1+(NBEL-NNN)/LRV
  94.       ENDIF
  95.       KPACKD=1
  96.       KPACKF=NPACK
  97. *
  98. *     boucle sur les packets, on active la structure petrov
  99. *
  100.       SEGINI PETROV
  101.       DO 7001 KPACK=KPACKD,KPACKF
  102. C
  103. C ======= A L'INTERIEUR DE CHAQUE PAQUET DE LRV ELEMENTS =======
  104. C
  105. C 1. Calcul des limites du paquet courant.
  106.          KDEB=1+(KPACK-1)*LRV
  107.          KFIN=MIN(NBEL,KDEB+LRV-1)
  108. C
  109.          DO 7002 K=KDEB,KFIN
  110.             KP=K-KDEB+1
  111.             NK = K + K0
  112.             NKR=(1-IK1)*(NK-1)+1
  113.             COEFK(KP)=IZTGG1.VPOCHA(NKR,1) + 1.D-30
  114. *            WRITE(*,*) IZTGG1.VPOCHA(NKR,1), COEFK(KP)
  115.             ANUK(KP)=ANUKK/COEFK(KP)
  116.  7002    CONTINUE
  117. C
  118. C        Une seule composante de projection car div se projete
  119. C        sur un scalaire.
  120. C
  121. c         write (*,*) 'IPADU', IPADU.LECT(/1)
  122.  
  123. C        CB215821 : DT n'est pas utilise mais doit etre initialise sinon NAN
  124.          DT=0.D0
  125.  
  126.          CALL SUPGEF(FN,GR,PG,XYZ,HR,PGSQ,RPG,AJ,
  127.      &   NES,NP,NPG,IAXI,XCOOR,
  128.      &   MELEME.NUM,KDEB,KFIN,LRV,IDCEN,CMD,IKOMP,
  129.      &   TETAN.VPOCHA,IPADU.LECT,IZTU1.VPOCHA,IPADU.LECT,NPTU,ANUK,
  130.      &   WT(1,1,1,1),WS(1,1,1,1),HK,PGSK,RPGK,AJK,AIRE,UIL,DUIL,
  131.      &   DTM1,DT,DTT1,DTT2,DIAEL,NUEL)
  132.  
  133.       DO 30 KE=KDEB,KFIN
  134.          NK=KE-KDEB+1
  135.  
  136. *      DO 12 I=1,NP
  137. *      J=NUM(I,KE)
  138. *      DO 12 N=1,IDIM
  139. *      XYZ(N,I)=XCOOR((J-1)*(IDIM+1)    +N)
  140. * 12   CONTINUE
  141.  
  142. *      CALL CALJBR
  143. *     &(FN,GR,PG,XYZ,HR,PGSQ,RPG,NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
  144. C     write(6,*)' Retour caljbr ',mp1,np,npg
  145.  
  146.  
  147.  
  148.          DO 324 K=1,IDIM
  149.             IPM4=KIPM(K)
  150.             DO     M=1,MP1
  151.                DO     I=1,NP
  152.                   U=0.D0
  153.                   DO 333 L=1,NPG
  154. *                    ON modifie la formule GBM 18/10/99
  155.              U=U+WT(NK,M,L,1)*HK(NK,K,I,L)*PGSK(NK,L)*DEUPI*RPGK(NK,L)
  156. *            U=U+FN(M,L)*HR(K,I,L)*PGSQ(L)*DEUPI*RPG(L)
  157.  333              CONTINUE
  158.  
  159.                   IF(IAXI.NE.0.AND.K.EQ.1)THEN
  160.                      DO 334 L=1,NPG
  161. *                      modifié par GBM 18/10/99
  162.                          U=U+WT(NK,M,L,1)*FN(I,L)*PGSK(NK,L)*DEUPI
  163. *                        U=U+FN(M,L)*FN(I,L)*PGSQ(L)*DEUPI
  164.  334                 CONTINUE
  165.                   ENDIF
  166. C
  167. C                 K1=1+(1-IK1)*(NK-1)
  168. C                 attention valable que pour coef au centre
  169.                   U=U*COEFK(NK)
  170. *                  WRITE(*,*) NK,M,I,COEFK(NK)
  171. C
  172.                   if(ikas.ne.2)then
  173.                      IPM4.AM(KE,I,M)=IPM4.AM(KE,I,M)+U
  174.                   else
  175.                      IPM4.AM(KE,M,I)=IPM4.AM(KE,M,I)+U
  176.                   endif
  177.  
  178. *            write(*,*) KE,PGSQ
  179.             ENDDO
  180.             ENDDO
  181.  324     CONTINUE
  182.  30   CONTINUE
  183.  7001 CONTINUE
  184.  
  185.       SEGSUP IZHR,IZFFM,PETROV
  186.  
  187.  
  188. C CAS MACRO CENTRE
  189.  
  190.       ELSEIF(KPRE.EQ.2)THEN
  191. * initialisation de npg au pif par TC
  192.         NPG=1
  193.         SEGINI PETROV
  194.         NOM0=NOMS(ITYPEL)//'    '
  195.         CALL KALPBG(NOM0,'FONFORM ',IZFFM)
  196.  
  197.         SEGACT IZFFM*MOD
  198.         IZHR=KZHR(1)
  199.         SEGACT IZHR*MOD
  200.         NES=GR(/1)
  201.         NPG=GR(/3)
  202.         IZF1=KTP(1)
  203.         SEGACT IZF1*MOD
  204.         MPG=IZF1.FN(/2)
  205.         NP=GR(/2)
  206.  
  207.         DO 40 KE=1,NEL
  208.  
  209.         NK=K0+KE
  210.  
  211.         IX=0
  212.         DO    I=1,NP
  213.         J=MELEME.NUM(I,KE)
  214.         DO    N=1,IDIM
  215.         IX=IX+1
  216.         XYZ1(IX)=XCOOR((J-1)*(IDIM+1)    +N)
  217.         ENDDO
  218.         ENDDO
  219.  
  220.       CALL CALJBR(FN,GR,PG,XYZ1,HR,PGSQ,RPG,NES,
  221.      & IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
  222.  
  223.       DO 424 K=1,IDIM
  224.       IPM4=KIPM(K)
  225.  
  226.       DO 423 I=1,NP
  227.       U=0.D0
  228.       DO 433 L=1,NPG
  229.       U=U+HR(K,I,L)*PGSQ(L)*DEUPI*RPG(L)
  230.  433  CONTINUE
  231.  
  232.       IF(IAXI.NE.0.AND.K.EQ.1)THEN
  233.       DO 434 L=1,NPG
  234.       U=U+FN(I,L)*PGSQ(L)*DEUPI
  235.  434  CONTINUE
  236.       ENDIF
  237.  
  238.       K1=1+(1-IK1)*(NK-1)
  239.       U=U*IZTGG1.VPOCHA(K1,1)
  240.  
  241.       if(ikas.ne.2)then
  242.       IPM4.AM(KE,I,1)=IPM4.AM(KE,I,1)+U
  243.       else
  244.       IPM4.AM(KE,1,I)=IPM4.AM(KE,1,I)+U
  245.       endif
  246.  
  247.  423  CONTINUE
  248.  424  CONTINUE
  249.  
  250.  40   CONTINUE
  251.  
  252.       SEGSUP IZHR,IZFFM,PETROV
  253.  
  254.       ENDIF
  255.  
  256.       RETURN
  257.  1002 FORMAT(10(1X,1PE11.4))
  258.  1040 FORMAT(1X,'CALCUL MATRICE AM ',I4/10(1X,1PE11.4))
  259.       END
  260.  
  261.  
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268.  
  269.  
  270.  
  271.  
  272.  
  273.  
  274.  
  275.  
  276.  
  277.  
  278.  
  279.  
  280.  
  281.  
  282.  
  283.  
  284.  
  285.  
  286.  

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