Télécharger la plaquette Cast3M


fr en

Télécharger xlap1d.eso

Retour à la liste

Numérotation des lignes : 

xlap1d

  1. C XLAP1D    SOURCE    CB215821  20/11/25    13:43:14     10792          
  2.       SUBROUTINE XLAP1D(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         : XLAP1D
  15. C DESCRIPTION : Calcul de la matrice jacobienne du résidu du laplacien
  16. C               VF 2D.
  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 e_t )
  24. C               Les contributions sont plus "compliquées" à calculer que
  25. C               les simples (cf. xlap1c) car on a à dériver des produits
  26. C               de fonctions de la vitesse.
  27. C               xlap1e 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       : XLAP1A : Calcul de la matrice jacobienne du
  45. C                    résidu du laplacien VF 2D.
  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 2D.
  80. C SORTIES            : -
  81. C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
  82. C***********************************************************************
  83. C VERSION    : v1, 10/12/2001, version initiale
  84. C HISTORIQUE : v1, 10/12/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 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
  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. *
  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,NLCGAU,NLCDRO
  151. *
  152.       REAL*8 ALPHA,UMALPH
  153.       REAL*8 DRD,DRG
  154.       REAL*8 DUXXF,DUXYF,DUYXF,DUYYF
  155.       REAL*8 XD,XF,XG,XFMXD,XFMXG
  156.       REAL*8 YD,YF,YG,YFMYD,YFMYG
  157.       REAL*8 ALPHAX,ALPHAY,CNX,CNY
  158.       REAL*8 SIGNOR,SURFFA,VOLUEL
  159.       REAL*8 RHOP,UP,VP
  160.       REAL*8 FACTOR,FACTDU,FACTDV
  161.       REAL*8 DUDRHO,DUDROU,DVDRHO,DVDROV
  162.       REAL*8 DSTDU
  163.       REAL*8 TAUXX,TAUXY,TAUYY
  164.       REAL*8 EPSIXX,EPSIXY,EPSIYX,EPSIYY
  165.       REAL*8 GAMMXD,GAMMXG,GAMMYD,GAMMYG
  166. *
  167. * Executable statements
  168. *
  169.       IF (IMPR.GT.2) WRITE(IOIMP,*) 'Entrée dans xlap1d.eso'
  170. *   On calcule les contributions à (d Res_{\rho e_t} / d var) ; var
  171. *   prenant successivement les valeurs :
  172. *      \rho,   \rho u,    \rho v,     \rho e_t
  173. *   On dérive les termes :  (\tens{\tau(\grad{u})} \prod \vect{u})
  174. *                            \pscal \vect{n}
  175. *   ce qui donne deux contributions.
  176. *   Contribution 1 :
  177. *     ( (d \tens{\tau} / d var) \prod \vect{u}) \pscal \vect{n}
  178. *   Contribution 2 :
  179. *     ( \tens{\tau}  \prod (d \vect{u} / d var)) \pscal \vect{n}
  180. *   Note :
  181. *     pas de contribution à (d Res_{\rho e_t} / d \rho e_t).
  182. *   Les noms de matrices élémentaires (type MATSIM) associées sont :
  183. *   RETRHO, RETROU, RETROV
  184.       IF (LCLIMV) THEN
  185.          SEGACT KRVIMP
  186.       ENDIF
  187.       IF (LCLITO) THEN
  188.          SEGACT KRTOIM
  189.          SEGACT MPTOIM
  190.       ENDIF
  191.       SEGACT NOMINC
  192.       SEGACT KRCENT
  193.       SEGACT KRFACE
  194.       SEGACT MELEFL
  195.       SEGACT MPSURF
  196.       SEGACT MPNORM
  197.       SEGACT MPVOLU
  198.       SEGACT MPMUC
  199.       SEGACT MPROC
  200.       SEGACT MPVITC
  201.       SEGACT MPGRVF
  202.       SEGACT ICOGRV
  203.       NSOUCH=ICOGRV.IMACHE(/1)
  204.       DO 1 ISOUCH=1,NSOUCH
  205.          MCOGRV=ICOGRV.IMACHE(ISOUCH)
  206.          JCOGRV=ICOGRV.ICHAML(ISOUCH)
  207.          SEGACT JCOGRV
  208.          KDUNDX=JCOGRV.IELVAL(1)
  209.          KDUNDY=JCOGRV.IELVAL(2)
  210.          SEGDES JCOGRV
  211.          SEGACT KDUNDX
  212.          SEGACT KDUNDY
  213.          SEGACT MCOGRV
  214.          NELEM=MCOGRV.NUM(/2)
  215.          NPTEL=MCOGRV.NUM(/1)
  216.          NPP1=NPTEL
  217.          NPP2=NPTEL+1
  218.          NEL1=NELEM
  219.          NEL2=NELEM
  220.          IEL1=1
  221.          IEL2=1
  222.          SEGINI RETRHO
  223.          SEGINI RETROU
  224.          SEGINI RETROV
  225.          SEGINI CCGRV2
  226.          RETRHO.NOMPRI(1:4)=NOMINC.MOTS(1)
  227.          RETRHO.NOMPRI(5:8)='    '
  228.          RETRHO.NOMDUA(1:4)=NOMINC.MOTS(4)
  229.          RETRHO.NOMDUA(5:8)='    '
  230.          RETROU.NOMPRI(1:4)=NOMINC.MOTS(2)
  231.          RETROU.NOMPRI(5:8)='    '
  232.          RETROU.NOMDUA(1:4)=NOMINC.MOTS(4)
  233.          RETROU.NOMDUA(5:8)='    '
  234.          RETROV.NOMPRI(1:4)=NOMINC.MOTS(3)
  235.          RETROV.NOMPRI(5:8)='    '
  236.          RETROV.NOMDUA(1:4)=NOMINC.MOTS(4)
  237.          RETROV.NOMDUA(5:8)='    '
  238.          DO 12 IELEM=1,NELEM
  239. *     Le premier point du support de ICOGRV est un point FACE
  240.             NGFACE=MCOGRV.NUM(1,IELEM)
  241.             NLFACE=KRFACE.LECT(NGFACE)
  242.             IF (NLFACE.EQ.0) THEN
  243.                WRITE(IOIMP,*) 'Erreur de programmation n°1'
  244.                GOTO 9999
  245.             ENDIF
  246. *     On calcule la contribution 2 à la matrice jacobienne IJACO de la
  247. *     face NGFACE
  248. *           (points duaux     : centres à gauche et à droite de la face)
  249. *           (points primaux   : une partie (bicoz conditions aux limites)
  250. *            de ceux du stencil pour le calcul du gradient
  251. *            à la face, ils doivent être des points centres
  252. *            ET les centres à gauche et à droite qui servent pour le
  253. *            calcul de la vitesse sur la face)
  254. *       Si la vitesse sur la face est imposée par les
  255. *       conditions aux limites, la contribution 2 de la face à IJACO est
  256. *       nulle.
  257.             LCTR2A=.TRUE.
  258.             IF (LCLIMV) THEN
  259.                NLCLV=KRVIMP.LECT(NGFACE)
  260.                IF (NLCLV.NE.0) THEN
  261.                   LCTR2A=.FALSE.
  262.                ENDIF
  263.             ENDIF
  264.             IF (LCTR2A) THEN
  265.                NGCGAU=MELEFL.NUM(1,NLFACE)
  266.                NGCDRO=MELEFL.NUM(3,NLFACE)
  267.                NLCGAU=KRCENT.LECT(NGCGAU)
  268.                NLCDRO=KRCENT.LECT(NGCDRO)
  269.                LMUR=(NGCGAU.EQ.NGCDRO)
  270. *     On calcule ensuite les valeurs des coefficients qui servent pour
  271. *     interpoler la vitesse sur la face NGFACE à partir des vitesses
  272. *     gauche (et droite si non mur) et des dérivées de la vitesse sur la
  273. *     face NGFACE (de la même manière que dans xlap12.eso)
  274.                IF (LMUR) THEN
  275. *       Parametres geometriques
  276.                   ICOORX=((IDIM+1)*(NGFACE-1))+1
  277.                   XF=MCOORD.XCOOR(ICOORX)
  278.                   YF=MCOORD.XCOOR(ICOORX+1)
  279.                   ICOORX=((IDIM+1)*(NGCGAU-1))+1
  280.                   XG=MCOORD.XCOOR(ICOORX)
  281.                   YG=MCOORD.XCOOR(ICOORX+1)
  282.                   XFMXG=XF-XG
  283.                   YFMYG=YF-YG
  284.                   ALPHA=0.0D0
  285.                   UMALPH=1.0D0
  286.                   GAMMXG=1.D0
  287.                   EPSIXX=XFMXG
  288.                   EPSIXY=YFMYG
  289.                   GAMMYG=1.D0
  290.                   EPSIYX=XFMXG
  291.                   EPSIYY=YFMYG
  292.                ELSEIF (.NOT.LMUR) THEN
  293.                   ICOORX=((IDIM+1)*(NGFACE-1))+1
  294.                   XF=MCOORD.XCOOR(ICOORX)
  295.                   YF=MCOORD.XCOOR(ICOORX+1)
  296.                   ICOORX=((IDIM+1)*(NGCGAU-1))+1
  297.                   XG=MCOORD.XCOOR(ICOORX)
  298.                   YG=MCOORD.XCOOR(ICOORX+1)
  299.                   XFMXG=XF-XG
  300.                   YFMYG=YF-YG
  301.                   DRG=SQRT((XFMXG*XFMXG)+(YFMYG*YFMYG))
  302.                   ICOORX=((IDIM+1)*(NGCDRO-1))+1
  303.                   XD=MCOORD.XCOOR(ICOORX)
  304.                   YD=MCOORD.XCOOR(ICOORX+1)
  305.                   XFMXD=XF-XD
  306.                   YFMYD=YF-YD
  307.                   DRD=SQRT((XFMXD*XFMXD)+(YFMYD*YFMYD))
  308.                   ALPHA=DRG/(DRG+DRD)
  309.                   UMALPH=1.0D0-ALPHA
  310.                   GAMMXG=UMALPH
  311.                   GAMMXD=ALPHA
  312.                   EPSIXX=(XFMXG*UMALPH)+(XFMXD*ALPHA)
  313.                   EPSIXY=(YFMYG*UMALPH)+(YFMYD*ALPHA)
  314.                   GAMMYG=UMALPH
  315.                   GAMMYD=ALPHA
  316.                   EPSIYX=(XFMXG*UMALPH)+(XFMXD*ALPHA)
  317.                   EPSIYY=(YFMYG*UMALPH)+(YFMYD*ALPHA)
  318.                ELSE
  319.                   WRITE(IOIMP,*) 'Erreur de programmation n°2'
  320.                   GOTO 9999
  321.                ENDIF
  322. *     On distingue le cas où la face est un bord du maillage (mur)
  323. *                  du cas où la face est interne au maillage
  324.                IF (.NOT.LMUR) THEN
  325.                   NPD=2
  326.                   NPP=NPTEL+1
  327.                ELSE
  328.                   NPD=1
  329.                   NPP=NPTEL
  330.                ENDIF
  331.                MU=UMALPH*MPMUC.VPOCHA(NLCGAU,1) +
  332.      &              ALPHA*MPMUC.VPOCHA(NLCDRO,1)
  333. *     On calcule tout d'abord les valeurs du tenseur des contraintes sur
  334. *     la face NGFACE de la même manière que dans xlap12.eso
  335.                IF (LCLITO) THEN
  336.                   NLFTOI=KRTOIM.LECT(NGFACE)
  337.                ELSE
  338.                   NLFTOI=0
  339.                ENDIF
  340.                IF (NLFTOI.NE.0) THEN
  341.                   TAUXX=MPTOIM.VPOCHA(NLFTOI,1)
  342.                   TAUYY=MPTOIM.VPOCHA(NLFTOI,2)
  343.                   TAUXY=MPTOIM.VPOCHA(NLFTOI,3)
  344.                ELSE
  345.                   DUXXF=MPGRVF.VPOCHA(NLFACE,1)
  346.                   DUXYF=MPGRVF.VPOCHA(NLFACE,2)
  347.                   DUYXF=MPGRVF.VPOCHA(NLFACE,3)
  348.                   DUYYF=MPGRVF.VPOCHA(NLFACE,4)
  349.                   DSTDU=(2.0D0/3.0D0)*(DUXXF+DUYYF)
  350.                   TAUXX=MU*((2.0D0*DUXXF)-DSTDU)
  351.                   TAUXY=MU*(DUXYF+DUYXF)
  352.                   TAUYY=MU*((2.0D0*DUYYF)-DSTDU)
  353.                ENDIF
  354. *     IPD=1 : point à gauche du point NGFACE
  355. *     IPD=2 : point à droite du point NGFACE
  356.                DO 122 IPD=1,NPD
  357.                   NPDUAL=MELEFL.NUM((2*IPD)-1,NLFACE)
  358.                   IF (.NOT.LMUR) THEN
  359.                      CCGRV2.POIDU2(IPD,IEL2)=NPDUAL
  360.                   ELSE
  361.                      CCGRV2.POIDU1(IPD,IEL1)=NPDUAL
  362.                   ENDIF
  363.                   NLCEND=KRCENT.LECT(NPDUAL)
  364.                   IF (NLCEND.EQ.0) THEN
  365.                      WRITE(IOIMP,*) 'Erreur grave n°1'
  366.                      GOTO 9999
  367.                   ENDIF
  368.                   DO 124 IPP=1,NPP
  369.                      IF (IPP.EQ.NPTEL) THEN
  370.                         NPPRIM=NGCGAU
  371.                      ELSEIF (IPP.EQ.NPTEL+1) THEN
  372.                         NPPRIM=NGCDRO
  373.                      ELSEIF (IPP.GE.1.AND.IPP.LT.NPTEL) THEN
  374.                         NPPRIM=MCOGRV.NUM(IPP+1,IELEM)
  375.                      ELSE
  376.                         WRITE(IOIMP,*) 'Erreur grave n°2'
  377.                         GOTO 9999
  378.                      ENDIF
  379.                      LCTR2B=.TRUE.
  380.                      IF (LCLIMV) THEN
  381.                         NLCLV=KRVIMP.LECT(NPPRIM)
  382.                         IF (NLCLV.NE.0) THEN
  383.                            LCTR2B=.FALSE.
  384.                         ENDIF
  385.                      ENDIF
  386.                      IF (.NOT.LCTR2B) THEN
  387. *     Lorsque une contribution est nulle, on fixe artificiellement le
  388. *     point primal égal au point dual.
  389.                         IF (.NOT.LMUR) THEN
  390.                            CCGRV2.POIPR2(IPP,IEL2)=NPDUAL
  391.                            RETRHO.VALMA2(IPD,IPP,IEL2)=0.D0
  392.                            RETROU.VALMA2(IPD,IPP,IEL2)=0.D0
  393.                            RETROV.VALMA2(IPD,IPP,IEL2)=0.D0
  394.                         ELSE
  395.                            CCGRV2.POIPR1(IPP,IEL1)=NPDUAL
  396.                            RETRHO.VALMA1(IPD,IPP,IEL1)=0.D0
  397.                            RETROU.VALMA1(IPD,IPP,IEL1)=0.D0
  398.                            RETROV.VALMA1(IPD,IPP,IEL1)=0.D0
  399.                         ENDIF
  400.                      ELSE
  401. *     Les contributions 2 valent :
  402. *     (d Res_{\rho e_t})_d / (d var)_p =
  403. *     +/-1 (normale sortante, rentrante)  (1/V_d) * (S_f)
  404. *     * [ [ [  ((n_x * \tau_xx) + (n_y * \tau_yx))
  405. *            * ((\epsilon_xx * \alpha_x) + (\epsilon_xy * \alpha_y))
  406. *           ]
  407. *         * ((du)_p / (d var)_p)
  408. *         ]
  409. *       + [ [  ((n_x * \tau_xy) + (n_y * \tau_yy))
  410. *            * ((\epsilon_yx * \alpha_x) + (\epsilon_yy * \alpha_y))
  411. *           ]
  412. *         * ((dv)_p / (d var)_p)
  413. *         ]
  414. *       ]
  415. *     et (m étant le centre de gauche dans le cas mur
  416. *         et les centre gauche, puis droite dans le cas normal)
  417. *     (d Res_{\rho e_t})_d / (d var)_m =
  418. *     +/-1 (normale sortante, rentrante)  (1/V_d) * (S_f)
  419. *     * [ [ [  ((n_x * \tau_xx) + (n_y * \tau_yx))
  420. *            * (\gamma_x)
  421. *           ]
  422. *         * ((du)_p / (d var)_p)
  423. *         ]
  424. *       + [ [  ((n_x * \tau_xy) + (n_y * \tau_yy))
  425. *            * (\gamma_y)
  426. *           ]
  427. *         * ((dv)_p / (d var)_p)
  428. *         ]
  429. *       ]
  430. *     sachant que l'expression donnant l'interpolation de la
  431. *     vitesse (u,v) sur une face i est de la forme :
  432. *     u_i =   (\gamma_x,gauche * u_gauche)
  433. *           + (\gamma_x,droite * u_droite) (=0 dans le cas mur)
  434. *           + (\epsilon_xx,i * (du/dx)_i)
  435. *           + (\epsilon_xy,i * (du/dy)_i)
  436. *     v_i =   (\gamma_y,gauche * v_gauche)
  437. *           + (\gamma_y,droite * v_droite) (=0 dans le cas mur)
  438. *           + (\epsilon_yx,i * (dv/dx)_i)
  439. *           + (\epsilon_yy,i * (dv/dy)_i)
  440. *     (voir dans xlap12.eso et ci-dessus pour les valeurs des coeffs
  441. *      \gamma et \epsilon)
  442.                         NLCENP=KRCENT.LECT(NPPRIM)
  443.                         IF (NLCENP.EQ.0) THEN
  444.                            WRITE(IOIMP,*) 'Erreur grave n°3'
  445.                            GOTO 9999
  446.                         ENDIF
  447. *       normale sortante pour IPD=1, rentrante pour IPD=2
  448.                         SIGNOR=(-1.D0)**(IPD+1)
  449.                         VOLUEL=MPVOLU.VPOCHA(NLCEND,1)
  450.                         SURFFA=MPSURF.VPOCHA(NLFACE,1)
  451.                         CNX   =MPNORM.VPOCHA(NLFACE,1)
  452.                         CNY   =MPNORM.VPOCHA(NLFACE,2)
  453.                         RHOP  =MPROC.VPOCHA(NLCENP,1)
  454.                         UP    =MPVITC.VPOCHA(NLCENP,1)
  455.                         VP    =MPVITC.VPOCHA(NLCENP,2)
  456.                         FACTOR=SIGNOR*(1.D0/VOLUEL)*SURFFA
  457.                         IF (IPP.EQ.NPTEL) THEN
  458.                            FACTDU= ((CNX*TAUXX)+(CNY*TAUXY))
  459.      $                          * GAMMXG
  460.                            FACTDV= ((CNX*TAUXY)+(CNY*TAUYY))
  461.      $                          * GAMMYG
  462.                         ELSEIF (IPP.EQ.NPTEL+1) THEN
  463.                            FACTDU= ((CNX*TAUXX)+(CNY*TAUXY))
  464.      $                          * GAMMXD
  465.                            FACTDV= ((CNX*TAUXY)+(CNY*TAUYY))
  466.      $                          * GAMMYD
  467.                         ELSEIF (IPP.GE.1.AND.IPP.LT.NPTEL) THEN
  468.                            ALPHAX=KDUNDX.VELCHE(IPP+1,IELEM)
  469.                            ALPHAY=KDUNDY.VELCHE(IPP+1,IELEM)
  470.                            FACTDU= ((CNX*TAUXX)+(CNY*TAUXY))
  471.      $                          * ((EPSIXX*ALPHAX)+(EPSIXY*ALPHAY))
  472.                            FACTDV= ((CNX*TAUXY)+(CNY*TAUYY))
  473.      $                          * ((EPSIYX*ALPHAX)+(EPSIYY*ALPHAY))
  474.                         ELSE
  475.                            WRITE(IOIMP,*) 'Erreur grave n°3'
  476.                            GOTO 9999
  477.                         ENDIF
  478.                         DUDRHO=-UP /RHOP
  479.                         DUDROU=1.D0/RHOP
  480.                         DVDRHO=-VP /RHOP
  481.                         DVDROV=1.D0/RHOP
  482.                         IF (.NOT.LMUR) THEN
  483.                            CCGRV2.POIPR2(IPP,IEL2)=NPPRIM
  484.                            RETRHO.VALMA2(IPD,IPP,IEL2)=
  485.      $                          FACTOR*( (FACTDU*DUDRHO)
  486.      $                          +(FACTDV*DVDRHO))
  487.                            RETROU.VALMA2(IPD,IPP,IEL2)=
  488.      $                          FACTOR*  (FACTDU*DUDROU)
  489.                            RETROV.VALMA2(IPD,IPP,IEL2)=
  490.      $                          FACTOR*  (FACTDV*DVDROV)
  491.                         ELSE
  492.                            CCGRV2.POIPR1(IPP,IEL1)=NPPRIM
  493.                            RETRHO.VALMA1(IPD,IPP,IEL1)=
  494.      $                          FACTOR*( (FACTDU*DUDRHO)
  495.      $                          +(FACTDV*DVDRHO))
  496.                            RETROU.VALMA1(IPD,IPP,IEL1)=
  497.      $                          FACTOR*  (FACTDU*DUDROU)
  498.                            RETROV.VALMA1(IPD,IPP,IEL1)=
  499.      $                          FACTOR*  (FACTDV*DVDROV)
  500.                         ENDIF
  501.                      ENDIF
  502.  124              CONTINUE
  503.  122           CONTINUE
  504.                IF (.NOT.LMUR) THEN
  505.                   IEL2=IEL2+1
  506.                ELSE
  507.                   IEL1=IEL1+1
  508.                ENDIF
  509.             ENDIF
  510.  12      CONTINUE
  511.          NPP1=NPTEL
  512.          NPP2=NPTEL+1
  513.          NEL1=IEL1-1
  514.          NEL2=IEL2-1
  515.          SEGADJ RETRHO
  516.          SEGADJ RETROU
  517.          SEGADJ RETROV
  518.          SEGADJ CCGRV2
  519.          CCGRV2.LMATSI(**)=RETRHO
  520.          CCGRV2.LMATSI(**)=RETROU
  521.          CCGRV2.LMATSI(**)=RETROV
  522. * On accumule les matrices résultantes dans IJACO
  523.          CALL AJMTK(CCGRV2,IJACO,IMPR,IRET)
  524.          IF (IRET.NE.0) GOTO 9999
  525.          SEGSUP RETRHO
  526.          SEGSUP RETROU
  527.          SEGSUP RETROV
  528.          SEGSUP CCGRV2
  529. *
  530.          SEGDES MCOGRV
  531.          SEGDES KDUNDY
  532.          SEGDES KDUNDX
  533.  1    CONTINUE
  534.       SEGDES ICOGRV
  535.       SEGDES MPGRVF
  536.       SEGDES MPVITC
  537.       SEGDES MPMUC
  538.       SEGDES MPROC
  539.       SEGDES MPVOLU
  540.       SEGDES MPNORM
  541.       SEGDES MPSURF
  542.       SEGDES MELEFL
  543.       SEGDES KRFACE
  544.       SEGDES KRCENT
  545.       SEGDES NOMINC
  546.       IF (LCLITO) THEN
  547.          SEGDES KRTOIM
  548.          SEGDES MPTOIM
  549.       ENDIF
  550.       IF (LCLIMV) THEN
  551.          SEGDES KRVIMP
  552.       ENDIF
  553. *
  554. * Normal termination
  555. *
  556.       IRET=0
  557.       RETURN
  558. *
  559. * Format handling
  560. *
  561. *
  562. * Error handling
  563. *
  564.  9999 CONTINUE
  565.       IRET=1
  566.       WRITE(IOIMP,*) 'An error was detected in subroutine xlap1d'
  567.       RETURN
  568. *
  569. * End of subroutine XLAP1D
  570. *
  571.       END
  572.  
  573.  
  574.  
  575.  
  576.  
  577.  
  578.  
  579.  
  580.  
  581.  
  582.  
  583.  
  584.  

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