Télécharger la plaquette Cast3M


fr en

Télécharger xlap2d.eso

Retour à la liste

Numérotation des lignes : 

xlap2d

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

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