Télécharger la plaquette Cast3M


fr en

Télécharger xlap2c.eso

Retour à la liste

Numérotation des lignes : 

xlap2c

  1. C XLAP2C    SOURCE    CB215821  20/11/25    13:43:20     10792          
  2.       SUBROUTINE XLAP2C(ICOGRV,MPROC,MPVITC,
  3.      $     MPVOLU,MPNORM,MPSURF,MELEFL,
  4.      $     KRFACE,KRCENT,LCLIMV,KRVIMP,LCLITO,KRTOIM,
  5.      $     NOMINC,
  6.      $     MPMUC,
  7.      $     IJACO,
  8.      $     IMPR,IRET)
  9.       IMPLICIT INTEGER(I-N)
  10.       IMPLICIT REAL*8 (A-H,O-Z)
  11. C***********************************************************************
  12. C NOM         : XLAP2C
  13. C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
  14. C               VF 3D.
  15. C               Ici, on ne calcule que les contributions 'simples'
  16. C               à la matrice jacobienne faisant intervenir les
  17. C               coefficients pour le calcul des gradients de vitesse
  18. C               (ICOGRV).
  19. C               (contributions à (d Res_{\rho u} / d var) et (d
  20. C                Res_{\rho v} / d var) et (d Res_{\rho w} / d var)
  21. C                var prenant successivement les valeurs :
  22. C                \rho,   \rho u,    \rho v, \rho w    \rho e_t )
  23. C
  24. C
  25. C NOTE PERSO  : On pourrait programmer ça plus lisiblement en stockant
  26. C               les contributions dans un tableau de pointeurs (2
  27. C               indices, c'est possible ?) et en effectuant des produits
  28. C               matriciels (coeff. x matrices de dérivées).
  29. C               On n'y coupera pas si on veut regrouper 2D et 3D...
  30. C               On ne va pas le faire.
  31. C
  32. C
  33. C LANGAGE     : ESOPE
  34. C AUTEUR      : Stéphane GOUNAND (CEA/DEN/DM2S/SFME/LTMF)
  35. C               mél : gounand@semt2.smts.cea.fr
  36. C***********************************************************************
  37. C APPELES (UTIL)   : AJMTK : ajoute un objet de type MATSIM (non
  38. C                            standard) à un objet de type MATRIK.
  39. C APPELE PAR       : XLAP2A : Calcul de la matrice jacobienne du
  40. C                    résidu du laplacien VF 2D.
  41. C***********************************************************************
  42. C ENTREES            : ICOGRV (type MCHELM) : coefficients pour le
  43. C                        calcul du gradient de la vitesse aux
  44. C                        interfaces.
  45. C                      MPROC (type MPOVAL) : masse volumique par
  46. C                                            élément.
  47. C                      MPVITC (type MPOVAL) : vitesse par élément.
  48. C                      MPVOLU (type MPOVAL) : volume des éléments.
  49. C                      MPNORM (type MPOVAL) : normale aux faces.
  50. C                      MPSURF (type MPOVAL) : surface des faces.
  51. C                      MELEFL (type MELEME) : connectivités face-(centre
  52. C                                             gauche, centre droit).
  53. C                      KRFACE (type MLENTI) : tableau de repérage dans
  54. C                        le maillage des faces des éléments.
  55. C                      KRCENT (type MLENTI) : tableau de repérage dans
  56. C                        le maillage des centres des éléments.
  57. C                      LCLIMV (type logique) : .TRUE. => CL de Dirichlet
  58. C                        sur la vitesse.
  59. C                      KRVIMP (type MLENTI) : tableau de repérage dans
  60. C                        maillage des CL de Dirichlet sur la
  61. C                        vitesse.
  62. C                      LCLITO (type logique) : .TRUE. => CL de Dirichlet
  63. C                        sur le tenseur des contraintes.
  64. C                      KRTOIM (type MLENTI) : tableau de repérage dans
  65. C                        maillage des CL de Dirichlet sur le tenseur des
  66. C                        contraintes
  67. C                      NOMINC (type MLMOTS) : noms des inconnues.
  68. C                      MPMUC (type MPOVAL) : viscosité dynamique (SI).
  69. C ENTREES/SORTIES    : IJACO (type MATRIK) : matrice jacobienne du
  70. C                        résidu du laplacien VF 2D.
  71. C SORTIES            : -
  72. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  73. C***********************************************************************
  74. C VERSION    : v1, 05/03/2002, version initiale
  75. C HISTORIQUE : v1, 05/03/2002, création
  76. C HISTORIQUE :
  77. C HISTORIQUE :
  78. C***********************************************************************
  79. C Prière de PRENDRE LE TEMPS de compléter les commentaires
  80. C en cas de modification de ce sous-programme afin de faciliter
  81. C la maintenance !
  82. C***********************************************************************
  83.  
  84. -INC PPARAM
  85. -INC CCOPTIO
  86. -INC SMCOORD
  87. -INC SMCHPOI
  88.       POINTEUR MPMUC.MPOVAL
  89.       POINTEUR MPROC.MPOVAL ,MPVITC.MPOVAL
  90.       POINTEUR MPSURF.MPOVAL,MPNORM.MPOVAL,MPVOLU.MPOVAL
  91. -INC SMCHAML
  92.       POINTEUR ICOGRV.MCHELM,JCOGRV.MCHAML
  93.       POINTEUR KDUNDX.MELVAL,KDUNDY.MELVAL,KDUNDZ.MELVAL
  94. -INC SMELEME
  95.       POINTEUR MELEFL.MELEME
  96.       POINTEUR MCOGRV.MELEME
  97. -INC SMLENTI
  98.       POINTEUR KRVIMP.MLENTI,KRTOIM.MLENTI
  99.       POINTEUR KRCENT.MLENTI,KRFACE.MLENTI
  100. -INC SMLMOTS
  101.       POINTEUR NOMINC.MLMOTS
  103.       POINTEUR IJACO.MATRIK
  104. *
  105. * Objet matrice élémentaire simplifié
  106. *
  107.       SEGMENT GMATSI
  108.       INTEGER POIPR1(NPP1,NEL1)
  109.       INTEGER POIDU1(1,NEL1)
  110.       INTEGER POIPR2(NPP2,NEL2)
  111.       INTEGER POIDU2(2,NEL2)
  112.       POINTEUR LMATSI(0).MATSIM
  113.       ENDSEGMENT
  114. *  Contributions simples de la part du gradient de
  115. *  vitesse (CTSGRV)
  116.       POINTEUR CTSGRV.GMATSI
  117.       SEGMENT MATSIM
  118.       CHARACTER*8 NOMPRI,NOMDUA
  119.       REAL*8 VALMA1(1,NPP1,NEL1)
  120.       REAL*8 VALMA2(2,NPP2,NEL2)
  121.       ENDSEGMENT
  122.       POINTEUR ROURHO.MATSIM
  123.       POINTEUR ROVRHO.MATSIM
  124.       POINTEUR ROWRHO.MATSIM
  125.       POINTEUR ROUROU.MATSIM
  126.       POINTEUR ROVROU.MATSIM
  127.       POINTEUR ROWROU.MATSIM
  128.       POINTEUR ROUROV.MATSIM
  129.       POINTEUR ROVROV.MATSIM
  130.       POINTEUR ROWROV.MATSIM
  131.       POINTEUR ROUROW.MATSIM
  132.       POINTEUR ROVROW.MATSIM
  133.       POINTEUR ROWROW.MATSIM
  134. *
  135.       REAL*8 MU
  136. *
  137.       INTEGER IMPR,IRET
  138. *
  139.       LOGICAL LCLIMV,LCLITO
  140.       LOGICAL LMUR
  141.       LOGICAL LCTRB1,LCTRB2
  142. *
  143.       INTEGER IELEM,IPD,IPP,ISOUCH,IEL1,IEL2
  144.       INTEGER NELEM,NPD,NPP,NSOUCH,NEL1,NEL2,NPP1,NPP2
  145.       INTEGER NGCDRO,NGCGAU,NGFACE,NPPRIM,NPDUAL
  146.       INTEGER NLCENP,NLCEND,NLFACE,NLCLTO,NLCLV
  147.       INTEGER NPTEL
  148. *
  149.       REAL*8 ALPHAX,ALPHAY,ALPHAZ,CNX,CNY,CNZ
  150.       REAL*8 SIGNOR,SURFFA,VOLUEL
  151.       REAL*8 RHOP,UP,VP,WP
  152.       REAL*8 FACTOR
  153.       REAL*8 DUDRHO,DUDROU
  154.       REAL*8 DVDRHO,DVDROV
  155.       REAL*8 DWDRHO,DWDROW
  156.       REAL*8 BROUDU,BROUDV,BROUDW
  157.       REAL*8 BROVDU,BROVDV,BROVDW
  158.       REAL*8 BROWDU,BROWDV,BROWDW
  159. *
  160.       INTEGER ICOORX,NLCGAU,NLCDRO
  161.       REAL*8 XF,YF,ZF,XG,YG,ZG,XFMXG,YFMYG,ZFMZG,DRG
  162.      &     ,XD,YD,ZD,XFMXD,YFMYD,ZFMZD,DRD,ALPHA,UMALPH
  163. C
  164. *
  165. * Executable statements
  166. *
  167.       IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans xlap2c.eso'
  168. *   On calcule les contributions à (d Res_{\rho u} / d var), (d
  169. *   Res_{\rho v} / d var) et (d Res_{\rho w} / d var) ; var prenant
  170. *   successivement les valeurs :
  171. *      \rho,   \rho u,    \rho v,  \rho w,    \rho e_t
  172. *   Note :
  173. *     - (d Res_{\rho u} / d \rho e_t) = 0
  174. *     - (d Res_{\rho v} / d \rho e_t) = 0
  175. *     - (d Res_{\rho w} / d \rho e_t) = 0
  176. *   On dérive les termes : \tens{\tau} \prod \vect{n}
  177. *   Les noms de matrices élémentaires (type MATSIM) associées sont :
  178. *   ROURHO, ROUROU, ROUROV, ROUROW,
  179. *   ROVRHO, ROVROU, ROVROV, ROVROW,
  180. *   ROWRHO, ROWROU, ROWROV, ROWROW.
  181.       IF (LCLIMV) THEN
  182.          SEGACT KRVIMP
  183.       ENDIF
  184.       IF (LCLITO) THEN
  185.          SEGACT KRTOIM
  186.       ENDIF
  187.       SEGACT NOMINC
  188.       SEGACT KRCENT
  189.       SEGACT KRFACE
  190.       SEGACT MELEFL
  191.       SEGACT MPSURF
  192.       SEGACT MPNORM
  193.       SEGACT MPVOLU
  194.       SEGACT MPMUC
  195.       SEGACT MPROC
  196.       SEGACT MPVITC
  197.       SEGACT ICOGRV
  198.       NSOUCH=ICOGRV.IMACHE(/1)
  199.       DO 1 ISOUCH=1,NSOUCH
  200.          MCOGRV=ICOGRV.IMACHE(ISOUCH)
  201.          JCOGRV=ICOGRV.ICHAML(ISOUCH)
  202.          SEGACT JCOGRV
  203.          KDUNDX=JCOGRV.IELVAL(1)
  204.          KDUNDY=JCOGRV.IELVAL(2)
  205.          KDUNDZ=JCOGRV.IELVAL(3)
  206.          SEGDES JCOGRV
  207.          SEGACT KDUNDX
  208.          SEGACT KDUNDY
  209.          SEGACT KDUNDZ
  210.          SEGACT MCOGRV
  211.          NELEM=MCOGRV.NUM(/2)
  212.          NPTEL=MCOGRV.NUM(/1)
  213.          NPP1=NPTEL-1
  214.          NPP2=NPTEL-1
  215.          NEL1=NELEM
  216.          NEL2=NELEM
  217.          IEL1=1
  218.          IEL2=1
  219.          SEGINI ROURHO
  220.          SEGINI ROVRHO
  221.          SEGINI ROWRHO
  222.          SEGINI ROUROU
  223.          SEGINI ROVROU
  224.          SEGINI ROWROU
  225.          SEGINI ROUROV
  226.          SEGINI ROVROV
  227.          SEGINI ROWROV
  228.          SEGINI ROUROW
  229.          SEGINI ROVROW
  230.          SEGINI ROWROW
  231.          SEGINI CTSGRV
  232.          ROURHO.NOMPRI(1:4)=NOMINC.MOTS(1)
  233.          ROURHO.NOMPRI(5:8)='    '
  234.          ROURHO.NOMDUA(1:4)=NOMINC.MOTS(2)
  235.          ROURHO.NOMDUA(5:8)='    '
  236.          ROVRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
  237.          ROVRHO.NOMPRI(5:8)='    '
  238.          ROVRHO.NOMDUA(1:4)=NOMINC.MOTS(3)
  239.          ROVRHO.NOMDUA(5:8)='    '
  240.          ROWRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
  241.          ROWRHO.NOMPRI(5:8)='    '
  242.          ROWRHO.NOMDUA(1:4)=NOMINC.MOTS(4)
  243.          ROWRHO.NOMDUA(5:8)='    '
  244.          ROUROU.NOMPRI(1:4)=NOMINC.MOTS(2)
  245.          ROUROU.NOMPRI(5:8)='    '
  246.          ROUROU.NOMDUA(1:4)=NOMINC.MOTS(2)
  247.          ROUROU.NOMDUA(5:8)='    '
  248.          ROVROU.NOMPRI(1:4)=NOMINC.MOTS(2)
  249.          ROVROU.NOMPRI(5:8)='    '
  250.          ROVROU.NOMDUA(1:4)=NOMINC.MOTS(3)
  251.          ROVROU.NOMDUA(5:8)='    '
  252.          ROWROU.NOMPRI(1:4)=NOMINC.MOTS(2)
  253.          ROWROU.NOMPRI(5:8)='    '
  254.          ROWROU.NOMDUA(1:4)=NOMINC.MOTS(4)
  255.          ROWROU.NOMDUA(5:8)='    '
  256.          ROUROV.NOMPRI(1:4)=NOMINC.MOTS(3)
  257.          ROUROV.NOMPRI(5:8)='    '
  258.          ROUROV.NOMDUA(1:4)=NOMINC.MOTS(2)
  259.          ROUROV.NOMDUA(5:8)='    '
  260.          ROVROV.NOMPRI(1:4)=NOMINC.MOTS(3)
  261.          ROVROV.NOMPRI(5:8)='    '
  262.          ROVROV.NOMDUA(1:4)=NOMINC.MOTS(3)
  263.          ROVROV.NOMDUA(5:8)='    '
  264.          ROWROV.NOMPRI(1:4)=NOMINC.MOTS(3)
  265.          ROWROV.NOMPRI(5:8)='    '
  266.          ROWROV.NOMDUA(1:4)=NOMINC.MOTS(4)
  267.          ROWROV.NOMDUA(5:8)='    '
  268.          ROUROW.NOMPRI(1:4)=NOMINC.MOTS(4)
  269.          ROUROW.NOMPRI(5:8)='    '
  270.          ROUROW.NOMDUA(1:4)=NOMINC.MOTS(2)
  271.          ROUROW.NOMDUA(5:8)='    '
  272.          ROVROW.NOMPRI(1:4)=NOMINC.MOTS(4)
  273.          ROVROW.NOMPRI(5:8)='    '
  274.          ROVROW.NOMDUA(1:4)=NOMINC.MOTS(3)
  275.          ROVROW.NOMDUA(5:8)='    '
  276.          ROWROW.NOMPRI(1:4)=NOMINC.MOTS(4)
  277.          ROWROW.NOMPRI(5:8)='    '
  278.          ROWROW.NOMDUA(1:4)=NOMINC.MOTS(4)
  279.          ROWROW.NOMDUA(5:8)='    '
  280.          DO 12 IELEM=1,NELEM
  281. *     Le premier point du support de ICOGRV est un point FACE
  282.             NGFACE=MCOGRV.NUM(1,IELEM)
  283.             NLFACE=KRFACE.LECT(NGFACE)
  284.             IF (NLFACE.EQ.0) THEN
  285.                WRITE(IOIMP,*) 'Erreur de programmation n°1'
  286.                GOTO 9999
  287.             ENDIF
  288. *     On calcule la contribution à la matrice jacobienne IJACO de la face
  289. *     NGFAC (points duaux     : centres à gauche et à droite de la face)
  290. *           (points primaux   : une partie (bicoz conditions aux limites)
  291. *            de ceux du stencil pour le calcul du gradient
  292. *            à la face, ils doivent être des points centres)
  293. *       Si le tenseur des contraintes sur la face est imposé par les
  294. *       conditions aux limites, la contribution de la face à IJACO est
  295. *       nulle.
  296.             LCTRB1=.TRUE.
  297.             IF (LCLITO) THEN
  298.                NLCLTO=KRTOIM.LECT(NGFACE)
  299.                IF (NLCLTO.NE.0) THEN
  300.                   LCTRB1=.FALSE.
  301.                ENDIF
  302.             ENDIF
  303.             IF (LCTRB1) THEN
  304.                NGCGAU=MELEFL.NUM(1,NLFACE)
  305.                NGCDRO=MELEFL.NUM(3,NLFACE)
  306.                NLCGAU=KRCENT.LECT(NGCGAU)
  307.                NLCDRO=KRCENT.LECT(NGCDRO)
  308.                LMUR=(NGCGAU.EQ.NGCDRO)
  309. *     On distingue le cas où la face est un bord du maillage (mur)
  310. *                  du cas où la face est interne au maillage
  311.                IF (.NOT.LMUR) THEN
  312.                   NPD=2
  313.                   ICOORX = ((IDIM + 1) * (NGFACE - 1))+1
  314.                   XF = MCOORD.XCOOR(ICOORX)
  315.                   YF = MCOORD.XCOOR(ICOORX+1)
  316.                   ZF = MCOORD.XCOOR(ICOORX+2)
  317.                   ICOORX = ((IDIM + 1) * (NGCGAU - 1))+1
  318.                   XG = MCOORD.XCOOR(ICOORX)
  319.                   YG = MCOORD.XCOOR(ICOORX+1)
  320.                   ZG = MCOORD.XCOOR(ICOORX+2)
  321.                   XFMXG = XF - XG
  322.                   YFMYG = YF - YG
  323.                   ZFMZG = ZF - ZG
  324.                   DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG)+(ZFMZG*ZFMZG))
  325.                   ICOORX = ((IDIM + 1) * (NGCDRO - 1))+1
  326.                   XD = MCOORD.XCOOR(ICOORX)
  327.                   YD = MCOORD.XCOOR(ICOORX+1)
  328.                   ZD = MCOORD.XCOOR(ICOORX+2)
  329.                   XFMXD = XF - XD
  330.                   YFMYD = YF - YD
  331.                   ZFMZD = ZF - ZD
  332.                   DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD)+(ZFMZD*ZFMZD))
  333.                   ALPHA=DRG/(DRG+DRD)
  334.                   UMALPH= 1.0D0 - ALPHA
  335.                ELSE
  336.                   NPD=1
  337.                   ALPHA=0.0D0
  338.                   UMALPH=1.0D0
  339.                ENDIF
  340.                MU=UMALPH*MPMUC.VPOCHA(NLCGAU,1) +
  341.      &              ALPHA*MPMUC.VPOCHA(NLCDRO,1)
  342.                NPP=NPTEL-1
  343. *     IPD=1 : point à gauche du point NGFACE
  344. *     IPD=2 : point à droite du point NGFACE
  345.                DO 122 IPD=1,NPD
  346.                   NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
  347.                   IF (.NOT.LMUR) THEN
  348.                      CTSGRV.POIDU2(IPD,IEL2)=NPDUAL
  349.                   ELSE
  350.                      CTSGRV.POIDU1(IPD,IEL1)=NPDUAL
  351.                   ENDIF
  352.                   NLCEND=KRCENT.LECT(NPDUAL)
  353.                   IF (NLCEND.EQ.0) THEN
  354.                      WRITE(IOIMP,*) 'Erreur grave n°1'
  355.                      GOTO 9999
  356.                   ENDIF
  357.                   DO 124 IPP=1,NPP
  358.                      NPPRIM=MCOGRV.NUM(IPP+1,IELEM)
  359.                      LCTRB2=.TRUE.
  360.                      IF (LCLIMV) THEN
  361.                         NLCLV=KRVIMP.LECT(NPPRIM)
  362.                         IF (NLCLV.NE.0) THEN
  363.                            LCTRB2=.FALSE.
  364.                         ENDIF
  365.                      ENDIF
  366.                      IF (.NOT.LCTRB2) THEN
  367. *     Lorsque une contribution est nulle, on fixe artificiellement le
  368. *     point primal égal au point dual.
  369.                         IF (.NOT.LMUR) THEN
  370.                            CTSGRV.POIPR2(IPP,IEL2)=NPDUAL
  371.                            ROURHO.VALMA2(IPD,IPP,IEL2)=0.D0
  372.                            ROVRHO.VALMA2(IPD,IPP,IEL2)=0.D0
  373.                            ROWRHO.VALMA2(IPD,IPP,IEL2)=0.D0
  374.                            ROUROU.VALMA2(IPD,IPP,IEL2)=0.D0
  375.                            ROVROU.VALMA2(IPD,IPP,IEL2)=0.D0
  376.                            ROWROU.VALMA2(IPD,IPP,IEL2)=0.D0
  377.                            ROUROV.VALMA2(IPD,IPP,IEL2)=0.D0
  378.                            ROVROV.VALMA2(IPD,IPP,IEL2)=0.D0
  379.                            ROWROV.VALMA2(IPD,IPP,IEL2)=0.D0
  380.                            ROUROW.VALMA2(IPD,IPP,IEL2)=0.D0
  381.                            ROVROW.VALMA2(IPD,IPP,IEL2)=0.D0
  382.                            ROWROW.VALMA2(IPD,IPP,IEL2)=0.D0
  383.                         ELSE
  384.                            CTSGRV.POIPR1(IPP,IEL1)=NPDUAL
  385.                            ROURHO.VALMA1(IPD,IPP,IEL1)=0.D0
  386.                            ROVRHO.VALMA1(IPD,IPP,IEL1)=0.D0
  387.                            ROWRHO.VALMA1(IPD,IPP,IEL1)=0.D0
  388.                            ROUROU.VALMA1(IPD,IPP,IEL1)=0.D0
  389.                            ROVROU.VALMA1(IPD,IPP,IEL1)=0.D0
  390.                            ROWROU.VALMA1(IPD,IPP,IEL1)=0.D0
  391.                            ROUROV.VALMA1(IPD,IPP,IEL1)=0.D0
  392.                            ROVROV.VALMA1(IPD,IPP,IEL1)=0.D0
  393.                            ROWROV.VALMA1(IPD,IPP,IEL1)=0.D0
  394.                            ROUROW.VALMA1(IPD,IPP,IEL1)=0.D0
  395.                            ROVROW.VALMA1(IPD,IPP,IEL1)=0.D0
  396.                            ROWROW.VALMA1(IPD,IPP,IEL1)=0.D0
  397.                         ENDIF
  398.                      ELSE
  399. *     Les contributions valent :
  400. * 1.  (d Res_{\rho u})_d / (d var)_p =
  401. *     +/-1 (normale sortante, rentrante)  (1/V_d) * (S_f) * \mu
  402. *     * [ [ ( 4/3 (n_x * \alpha_x) + (n_y * \alpha_y)
  403. *                                  + (n_z * \alpha_z)) *
  404. *           ((du)_p / (d var)_p)]
  405. *       + [ (-2/3 (n_x * \alpha_y) + (n_y * \alpha_x)) *
  406. *           ((dv)_p / (d var)_p)]
  407. *       + [ (-2/3 (n_x * \alpha_z) + (n_z * \alpha_x)) *
  408. *           ((dw)_p / (d var)_p)]
  409. *       ]
  410. * 2.  (d Res_{\rho v})_d / (d var)_p =
  411. *     +/-1 (normale sortante, rentrante)  (1/V_d) * (S_f) * \mu
  412. *     * [ [ (-2/3 (n_y * \alpha_x) + (n_x * \alpha_y)) *
  413. *           ((du)_p / (d var)_p)]
  414. *       + [ ( 4/3 (n_y * \alpha_y) + (n_x * \alpha_x)
  415. *                                  + (n_z * \alpha_z)) *
  416. *           ((dv)_p / (d var)_p)]
  417. *       + [ (-2/3 (n_y * \alpha_z) + (n_z * \alpha_y)) *
  418. *           ((dw)_p / (d var)_p)]
  419. *       ]
  420. * 3.  (d Res_{\rho w})_d / (d var)_p =
  421. *     +/-1 (normale sortante, rentrante)  (1/V_d) * (S_f) * \mu
  422. *     * [ [ (-2/3 (n_z * \alpha_x) + (n_x * \alpha_z)) *
  423. *           ((du)_p / (d var)_p)]
  424. *       + [ (-2/3 (n_z * \alpha_y) + (n_y * \alpha_z)) *
  425. *           ((dv)_p / (d var)_p)]
  426. *       + [ ( 4/3 (n_z * \alpha_z) + (n_x * \alpha_x)
  427. *                                  + (n_y * \alpha_y)) *
  428. *           ((dw)_p / (d var)_p)]
  429. *       ]
  430. *
  431. *     avec :
  432. *          (du)_p / (d \rho)_p   = - u / \rho_p
  433. *          (du)_p / (d \rho u)_p =   1 / \rho_p
  434. *          (du)_p / (d \rho v)_p =   0
  435. *          (du)_p / (d \rho w)_p =   0
  436. *          (dv)_p / (d \rho)_p   = - v / \rho_p
  437. *          (dv)_p / (d \rho u)_p =   0
  438. *          (dv)_p / (d \rho v)_p =   1 / \rho_p
  439. *          (dv)_p / (d \rho w)_p =   0
  440. *          (dv)_p / (d \rho)_p   = - w / \rho_p
  441. *          (dv)_p / (d \rho u)_p =   0
  442. *          (dv)_p / (d \rho v)_p =   0
  443. *          (dv)_p / (d \rho w)_p =   1 / \rho_p
  444. *
  445.                         NLCENP=KRCENT.LECT(NPPRIM)
  446.                         IF (NLCENP.EQ.0) THEN
  447.                            WRITE(IOIMP,*) 'Erreur grave n°2'
  448.                            GOTO 9999
  449.                         ENDIF
  450. *       normale sortante pour IPD=1, rentrante pour IPD=2
  451.                         SIGNOR=(-1.D0)**(IPD+1)
  452.                         VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
  453.                         SURFFA=MPSURF.VPOCHA(NLFACE,1)
  454.                         CNX   =MPNORM.VPOCHA(NLFACE,1)
  455.                         CNY   =MPNORM.VPOCHA(NLFACE,2)
  456.                         CNZ   =MPNORM.VPOCHA(NLFACE,3)
  457.                         ALPHAX=KDUNDX.VELCHE(IPP+1,IELEM)
  458.                         ALPHAY=KDUNDY.VELCHE(IPP+1,IELEM)
  459.                         ALPHAZ=KDUNDZ.VELCHE(IPP+1,IELEM)
  460.                         RHOP  =MPROC.VPOCHA(NLCENP,1)
  461.                         UP    =MPVITC.VPOCHA(NLCENP,1)
  462.                         VP    =MPVITC.VPOCHA(NLCENP,2)
  463.                         WP    =MPVITC.VPOCHA(NLCENP,3)
  464.                         FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA*MU
  465.                         BROUDU=(( 4.D0/3.D0)*(CNX*ALPHAX))
  466.      $                       + (CNY*ALPHAY) + (CNZ*ALPHAZ)
  467.                         BROUDV=((-2.D0/3.D0)*(CNX*ALPHAY))
  468.      $                       + (CNY*ALPHAX)
  469.                         BROUDW=((-2.D0/3.D0)*(CNX*ALPHAZ))
  470.      $                       + (CNZ*ALPHAX)
  471.                         BROVDU=((-2.D0/3.D0)*(CNY*ALPHAX))
  472.      $                       + (CNX*ALPHAY)
  473.                         BROVDV=(( 4.D0/3.D0)*(CNY*ALPHAY))
  474.      $                       + (CNX*ALPHAX) + (CNZ*ALPHAZ)
  475.                         BROVDW=((-2.D0/3.D0)*(CNY*ALPHAZ))
  476.      $                       + (CNZ*ALPHAY)
  477.                         BROWDU=((-2.D0/3.D0)*(CNZ*ALPHAX))
  478.      $                       + (CNX*ALPHAZ)
  479.                         BROWDV=((-2.D0/3.D0)*(CNZ*ALPHAY))
  480.      $                       + (CNY*ALPHAZ)
  481.                         BROWDW=(( 4.D0/3.D0)*(CNZ*ALPHAZ))
  482.      $                       + (CNX*ALPHAX) + (CNY*ALPHAY)
  483.                         DUDRHO=-UP /RHOP
  484.                         DUDROU=1.D0/RHOP
  485.                         DVDRHO=-VP /RHOP
  486.                         DVDROV=1.D0/RHOP
  487.                         DWDRHO=-WP /RHOP
  488.                         DWDROW=1.D0/RHOP
  489.                         IF (.NOT.LMUR) THEN
  490.                            CTSGRV.POIPR2(IPP,IEL2)=NPPRIM
  491.                            ROURHO.VALMA2(IPD,IPP,IEL2)=
  492.      $                          FACTOR*((BROUDU*DUDRHO)
  493.      $                          +(BROUDV*DVDRHO)+(BROUDW*DWDRHO))
  494.                            ROVRHO.VALMA2(IPD,IPP,IEL2)=
  495.      $                          FACTOR*((BROVDU*DUDRHO)
  496.      $                          +(BROVDV*DVDRHO)+(BROVDW*DWDRHO))
  497.                            ROWRHO.VALMA2(IPD,IPP,IEL2)=
  498.      $                          FACTOR*((BROWDU*DUDRHO)
  499.      $                          +(BROWDV*DVDRHO)+(BROWDW*DWDRHO))
  500.                            ROUROU.VALMA2(IPD,IPP,IEL2)=
  501.      $                          FACTOR* (BROUDU*DUDROU)
  502.                            ROVROU.VALMA2(IPD,IPP,IEL2)=
  503.      $                          FACTOR* (BROVDU*DUDROU)
  504.                            ROWROU.VALMA2(IPD,IPP,IEL2)=
  505.      $                          FACTOR* (BROWDU*DUDROU)
  506.                            ROUROV.VALMA2(IPD,IPP,IEL2)=
  507.      $                          FACTOR* (BROUDV*DVDROV)
  508.                            ROVROV.VALMA2(IPD,IPP,IEL2)=
  509.      $                          FACTOR* (BROVDV*DVDROV)
  510.                            ROWROV.VALMA2(IPD,IPP,IEL2)=
  511.      $                          FACTOR* (BROWDV*DVDROV)
  512.                            ROUROW.VALMA2(IPD,IPP,IEL2)=
  513.      $                          FACTOR* (BROUDW*DWDROW)
  514.                            ROVROW.VALMA2(IPD,IPP,IEL2)=
  515.      $                          FACTOR* (BROVDW*DWDROW)
  516.                            ROWROW.VALMA2(IPD,IPP,IEL2)=
  517.      $                          FACTOR* (BROWDW*DWDROW)
  518.                         ELSE
  519.                            CTSGRV.POIPR1(IPP,IEL1)=NPPRIM
  520.                            ROURHO.VALMA1(IPD,IPP,IEL1)=
  521.      $                          FACTOR*((BROUDU*DUDRHO)
  522.      $                          +(BROUDV*DVDRHO)+(BROUDW*DWDRHO))
  523.                            ROVRHO.VALMA1(IPD,IPP,IEL1)=
  524.      $                          FACTOR*((BROVDU*DUDRHO)
  525.      $                          +(BROVDV*DVDRHO)+(BROVDW*DWDRHO))
  526.                            ROWRHO.VALMA1(IPD,IPP,IEL1)=
  527.      $                          FACTOR*((BROWDU*DUDRHO)
  528.      $                          +(BROWDV*DVDRHO)+(BROWDW*DWDRHO))
  529.                            ROUROU.VALMA1(IPD,IPP,IEL1)=
  530.      $                          FACTOR* (BROUDU*DUDROU)
  531.                            ROVROU.VALMA1(IPD,IPP,IEL1)=
  532.      $                          FACTOR* (BROVDU*DUDROU)
  533.                            ROWROU.VALMA1(IPD,IPP,IEL1)=
  534.      $                          FACTOR* (BROWDU*DUDROU)
  535.                            ROUROV.VALMA1(IPD,IPP,IEL1)=
  536.      $                          FACTOR* (BROUDV*DVDROV)
  537.                            ROVROV.VALMA1(IPD,IPP,IEL1)=
  538.      $                          FACTOR* (BROVDV*DVDROV)
  539.                            ROWROV.VALMA1(IPD,IPP,IEL1)=
  540.      $                          FACTOR* (BROWDV*DVDROV)
  541.                            ROUROW.VALMA1(IPD,IPP,IEL1)=
  542.      $                          FACTOR* (BROUDW*DWDROW)
  543.                            ROVROW.VALMA1(IPD,IPP,IEL1)=
  544.      $                          FACTOR* (BROVDW*DWDROW)
  545.                            ROWROW.VALMA1(IPD,IPP,IEL1)=
  546.      $                          FACTOR* (BROWDW*DWDROW)
  547.                         ENDIF
  548.                      ENDIF
  549.  124              CONTINUE
  550.  122           CONTINUE
  551.                IF (.NOT.LMUR) THEN
  552.                   IEL2=IEL2+1
  553.                ELSE
  554.                   IEL1=IEL1+1
  555.                ENDIF
  556.             ENDIF
  557.  12      CONTINUE
  558.          NPP1=NPTEL-1
  559.          NPP2=NPTEL-1
  560.          NEL1=IEL1-1
  561.          NEL2=IEL2-1
  562.          SEGADJ ROURHO
  563.          SEGADJ ROVRHO
  564.          SEGADJ ROWRHO
  565.          SEGADJ ROUROU
  566.          SEGADJ ROVROU
  567.          SEGADJ ROWROU
  568.          SEGADJ ROUROV
  569.          SEGADJ ROVROV
  570.          SEGADJ ROWROV
  571.          SEGADJ ROUROW
  572.          SEGADJ ROVROW
  573.          SEGADJ ROWROW
  574.          SEGADJ CTSGRV
  575.          CTSGRV.LMATSI(**)=ROURHO
  576.          CTSGRV.LMATSI(**)=ROVRHO
  577.          CTSGRV.LMATSI(**)=ROWRHO
  578.          CTSGRV.LMATSI(**)=ROUROU
  579.          CTSGRV.LMATSI(**)=ROVROU
  580.          CTSGRV.LMATSI(**)=ROWROU
  581.          CTSGRV.LMATSI(**)=ROUROV
  582.          CTSGRV.LMATSI(**)=ROVROV
  583.          CTSGRV.LMATSI(**)=ROWROV
  584.          CTSGRV.LMATSI(**)=ROUROW
  585.          CTSGRV.LMATSI(**)=ROVROW
  586.          CTSGRV.LMATSI(**)=ROWROW
  587. * On accumule les matrices résultantes dans IJACO
  588.          CALL AJMTK(CTSGRV,IJACO,IMPR,IRET)
  589.          IF (IRET.NE.0) GOTO 9999
  590.          SEGSUP ROURHO
  591.          SEGSUP ROVRHO
  592.          SEGSUP ROWRHO
  593.          SEGSUP ROUROU
  594.          SEGSUP ROVROU
  595.          SEGSUP ROWROU
  596.          SEGSUP ROUROV
  597.          SEGSUP ROVROV
  598.          SEGSUP ROWROV
  599.          SEGSUP ROUROW
  600.          SEGSUP ROVROW
  601.          SEGSUP ROWROW
  602.          SEGSUP CTSGRV
  603. *
  604.          SEGDES MCOGRV
  605.          SEGDES KDUNDZ
  606.          SEGDES KDUNDY
  607.          SEGDES KDUNDX
  608.  1    CONTINUE
  609.       SEGDES ICOGRV
  610.       SEGDES MPVITC
  611.       SEGDES MPMUC
  612.       SEGDES MPROC
  613.       SEGDES MPVOLU
  614.       SEGDES MPNORM
  615.       SEGDES MPSURF
  616.       SEGDES MELEFL
  617.       SEGDES KRFACE
  618.       SEGDES KRCENT
  619.       SEGDES NOMINC
  620.       IF (LCLITO) THEN
  621.          SEGDES KRTOIM
  622.       ENDIF
  623.       IF (LCLIMV) THEN
  624.          SEGDES KRVIMP
  625.       ENDIF
  626. *
  627. * Normal termination
  628. *
  629.       IRET=0
  630.       RETURN
  631. *
  632. * Format handling
  633. *
  634. *
  635. * Error handling
  636. *
  637.  9999 CONTINUE
  638.       IRET=1
  639.       WRITE(IOIMP,*) 'An error was detected in subroutine xlap2c'
  640.       RETURN
  641. *
  642. * End of subroutine XLAP2C
  643. *
  644.       END
  645.  
  646.  
  647.  
  648.  
  649.  
  650.  
  651.  
  652.  
  653.  
  654.  
  655.  
  656.  
  657.  

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