Télécharger la plaquette Cast3M


fr en

Télécharger ylap2d.eso

Retour à la liste

Numérotation des lignes : 

ylap2d

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

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