Télécharger la plaquette Cast3M


fr en

Télécharger pente2.eso

Retour à la liste

Numérotation des lignes : 

pente2

  1. C PENTE2    SOURCE    CB215821  20/11/25    13:35:38     10792          
  2.       SUBROUTINE PENTE2(IOP2,INORM,MLECEN,MLEFAC,
  3.      &     MLENCL,IMCHAM,
  4.      &     NCOMP,MPOCHP,MPOVCL,MPOGRA,
  5.      &     MPOMIN,MPOMAX)
  6. C
  7. C************************************************************************
  8. C
  9. C PROJET            :  CASTEM 2000
  10. C
  11. C NOM               :  PENTE2
  12. C
  13. C DESCRIPTION       :  Cette subroutine est appellée par la subroutine
  14. C                      PENTE1 (calcul du gradient d'un CHPOINT de type
  15. C                      CENTRE)
  16. C                      Elle contient la partie du calcul du gradient.
  17. C
  18. C LANGUAGE          :  ESOPE 2000 avec extensions CISI
  19. C
  20. C AUTEUR            :  A. BECCANTINI, DRN/DMT/SEMT/TTMF
  21. C AUTEUR (modif.)   :  R. MOREL, DRN/DMT/SEMT/TTMF
  22. C
  23. C************************************************************************
  24. C
  25. C
  26. C APPELES (Outils)  : LICHT, ERREUR
  27. C
  28. C
  29. C************************************************************************
  30. C
  31. C ENTREES  : IOP2   :  Entier indiquant le type de reconstruction (cf PENT)
  32. C
  33. C            INORM   : Pointeur sur le CHPOINT des normales aux faces
  34. C
  35. C            MLECEN  : table numerotation global/local pour les
  36. C                            points CENTRE (segment MLENTI)
  37. C
  38. C            MLEFAC  : table numerotation global/local pour les
  39. C                            points FACE
  40. C
  41. C            MLENCL  : table numerotation global/local pour les
  42. C                            SPG du CHPOINT de C.L. ou
  43. C                      table de 0 si ce CHPOINT n'existe pas
  44. C
  45. C            IMCHAM  : MCHAML de type 'GRADGEO'
  46. C
  47. C            NCOMP   : Nombre de composantes du CHPOINT dont
  48. C                      on veut le gradient
  49. C
  50. C            MPOCHP  : MPOVAL du CHPOINT dont on veut le gradient
  51. C
  52. C            MPOVCL  : MPOVAL du CHPOINT de C.L.
  53. C
  54. C SORTIES  : MPOGRA  : MPOVAL du gradient
  55. C
  56. C            MPOMAX  : MPOVAL du max sur le stencil
  57. C
  58. C            MPOMIN  : MPOVAL du min sur le stencil
  59. C
  60. C
  61. C************************************************************************
  62. C
  63. C HISTORIQUE (Anomalies et modifications éventuelles)
  64. C
  65. C HISTORIQUE :  Cree le 4-6-1998
  66. C
  67. C HISTORIQUE :  Modifie le 20-10-1998 pour extension 3D
  68. C               Modifie le 28-04-2000 pour reconstruction quadatique
  69. C               exacte (A. Beccantini)
  70. C HISTORIQUE :  Modifie le 05.07.2003: elimination reconstruction
  71. C               quadratique exacte.
  72. C               Algorithme de BArth-Jespersen remplacé par une
  73. C               reconstruction lineaire exacte qui se base sur
  74. C               LSM
  75. C************************************************************************
  76. C
  77.       IMPLICIT INTEGER(I-N)
  78.       IMPLICIT REAL*8(A-H,O-Z)
  79.  
  80. -INC PPARAM
  81. -INC CCOPTIO
  82. -INC SMCHAML
  83. -INC SMELEME
  84. -INC SMCHPOI
  85. -INC SMLENTI
  86. C
  87.       POINTEUR INORM.MCHPOI
  88.       POINTEUR MLENCL.MLENTI, MLECEN.MLENTI, MLEFAC.MLENTI
  89.       POINTEUR MPOMAX.MPOVAL, MPOMIN.MPOVAL, MPONOR.MPOVAL,
  90.      &         MPOVCL.MPOVAL, MPOCHP.MPOVAL, MPOGRA.MPOVAL
  91.       POINTEUR MELVX.MELVAL,  MELVY.MELVAL,  MELVZ.MELVAL,
  92.      &         MELVXX.MELVAL, MELVYY.MELVAL, MELVZZ.MELVAL,
  93.      &         MELVXY.MELVAL, MELVXZ.MELVAL, MELVYZ.MELVAL
  94. C
  95.       INTEGER IOP2,NCOMP,NLELEM,NVOI,I3,IMCHAM
  96.      &     ,NGCE,NLCE,NTOTV,NGCV,NLCV,ICOM,IGEOM,K,NBELEM,NBNN
  97.      &     ,NLCF,NLCL,NMAI
  98.       REAL*8  VAL,VCOMP(3),VN
  99.       REAL*8  NX,NY,NZ
  100.       CHARACTER*8 TYPE
  101. C
  102. C**** Segments déjà activés
  103. C
  104. C     MPOVCL
  105. C     MLENCL
  106. C     MLECEN
  107. C     MLEFAC
  108. C     MPOCHP
  109. C     MPOGRA*MOD
  110. C     MPOMIN*MOD
  111. C     MPOMAX*MOD
  112. C
  113. C**** N.B. MPOMAX.VPOCHA, MPOMIN.VPOCHA déjà initialisé:
  114. C
  115. C     MPOMAX.VPOCHA = MPOCHP.VPOCHA
  116. C     MPOMIN.VPOCHA = MPOCHP.VPOCHA
  117. C
  118.       CALL LICHT(INORM,MPONOR,TYPE,IGEOM)
  119. C
  120. C**** En LICHT
  121. C     SEGACT*MOD MPONOR
  122. C     MPONOR: pointeur sur les normales aux faces
  123. C
  124.       MCHELM = IMCHAM
  125.       SEGACT MCHELM
  126.       NMAI=MCHELM.ICHAML(/1)
  127. C
  128. C*********************************************
  129. C Premier cas : on traite un champ scalaire en
  130. C               reconstruction non-quadratique
  131. C*********************************************
  132. C
  133. C
  134.       IF(IOP2.EQ.1) THEN
  135. C
  136. C**** N.B:  IOP2 = 1 -> EULESCAL
  137. C
  138.          DO 10 K=1,NMAI
  139. C
  140. C********** Boucle sur les differents maillages elementaires
  141. C           IPT1 = maillage duale
  142. C
  143.             IPT1=MCHELM.IMACHE(K)
  144.             SEGACT IPT1
  145.             NBNN=IPT1.NUM(/1)
  146.             NTOTV=NBNN
  147. C
  148. C********** NTOTV = Nombre de voisins
  149. C
  150.             NBELEM = IPT1.NUM(/2)
  151.             MCHAM1 = MCHELM.ICHAML(K)
  152.             SEGACT MCHAM1
  153.             MELVX = MCHAM1.IELVAL(1)
  154.             MELVY = MCHAM1.IELVAL(2)
  155.             SEGACT MELVX
  156.             SEGACT MELVY
  157.             IF (IDIM .EQ. 3) THEN
  158.                MELVZ = MCHAM1.IELVAL(3)
  159.                SEGACT MELVZ
  160.             ENDIF
  161.             SEGDES MCHAM1
  162.             DO NLELEM = 1, NBELEM, 1
  163. C
  164. C************* Boucle sur les ELTs
  165. C              NLELEM = Numero local d'ELT dans le maillage 1
  166. C
  167. C              NGCE = numero global centre d'elt
  168. C              NLCE = numero local centre d'elt
  169. C
  170. C              MLECEN  = table numerotation global/local pour les
  171. C                        points CENTRE (segment MLENTI)
  172. C
  173. C              MLEFAC  = table numerotation global/local pour les
  174. C                        points FACE
  175. C
  176. C              MLENCL  = table numerotation global/local pour les
  177. C                        SPG du CHPOINT de C.L. ou
  178. C                        table de 0 si ce CHPOINT n'existe pas
  179. C
  180.                NGCE = IPT1.NUM(NBNN,NLELEM)
  181.                NLCE = MLECEN.LECT(NGCE)
  182.                DO NVOI = 1, NTOTV
  183. C
  184. C**************** Boucle sur les voisins de NGCE
  185. C
  186. C                 NGCV  = Numero global voisin
  187. C                 NLCV  =   "    local   de NGCV dans le maillage CENTRE
  188. C                 NLCL  =   "    local   de NGCV dans le SPG du CHPOINT
  189. C                         des C.L.
  190. C
  191.                   NGCV = IPT1.NUM(NVOI,NLELEM)
  192.                   NLCV = MLECEN.LECT(NGCV)
  193.                   IF(NLCV.EQ.0)THEN
  194.                      NLCL = MLENCL.LECT(NGCV)
  195. C
  196. C**************** NGCV n'est pas un point centre, mais un point face,
  197. C                 et on est sur les bords
  198. C
  199.                      DO I3 = 1, NCOMP
  200. C
  201. C********************** I3   : numero de composante du CHPOI
  202. C                       ICOM : numero de composante du gradient
  203. C
  204.                         ICOM = IDIM*(I3 -1)+1
  205.                         IF (NLCL .EQ. 0) THEN
  206. C
  207. C************************* NGCV n'est pas un point du SPG de C.L.
  208. C
  209.                            VAL =  MPOCHP.VPOCHA(NLCE,I3)
  210.                         ELSE
  211. C
  212. C************************* NGCV est un point du SPG de C.L.
  213. C
  214.                            VAL = MPOVCL.VPOCHA(NLCL,I3)
  215.                         ENDIF
  216.                         MPOGRA.VPOCHA(NLCE,ICOM)= MPOGRA.VPOCHA(NLCE
  217.      $                       ,ICOM)+VAL * MELVX.VELCHE(NVOI,NLELEM)
  218.                         MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE
  219.      $                       ,ICOM+1) +VAL * MELVY.VELCHE(NVOI,NLELEM)
  220.                         IF (IDIM .EQ. 3) THEN
  221.                            MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE
  222.      $                          ,ICOM+2) +VAL * MELVZ.VELCHE(NVOI
  223.      $                          ,NLELEM)
  224.                         ENDIF
  225.                         MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3
  226.      $                       ),VAL)
  227.                         MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3
  228.      $                       ),VAL)
  229.                      ENDDO
  230.                   ELSE
  231. C
  232. C********** NGCV est un point centre
  233. C
  234.                      DO I3 = 1, NCOMP
  235.                         ICOM = IDIM*(I3 -1)+1
  236.                         VAL = MPOCHP.VPOCHA(NLCV,I3)
  237.                         MPOGRA.VPOCHA(NLCE,ICOM)=MPOGRA.VPOCHA(NLCE,ICOM
  238.      $                       )+VAL * MELVX.VELCHE(NVOI,NLELEM)
  239.                         MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE
  240.      $                       ,ICOM+1)+VAL * MELVY.VELCHE(NVOI,NLELEM)
  241.                         IF (IDIM .EQ. 3) THEN
  242.                            MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE
  243.      $                          ,ICOM+2)+VAL * MELVZ.VELCHE(NVOI,NLELEM
  244.      $                          )
  245.                         ENDIF
  246.                         MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3
  247.      $                       ),VAL)
  248.                         MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3
  249.      $                       ),VAL)
  250.                      ENDDO
  251.                   ENDIF
  252.                ENDDO
  253. C
  254. C******* Fin boucle sur les voisins de NGCE
  255. C
  256.             ENDDO
  257.             SEGDES MELVX
  258.             SEGDES MELVY
  259.             IF (IDIM .EQ. 3) SEGDES MELVZ
  260.             SEGDES IPT1
  261.  10      CONTINUE
  262. C
  263. C**********************************************************
  264. C Deuxieme cas : on traite un champ vectoriel (IOP2=2)
  265. C**********************************************************
  266. C
  267. C
  268.       ELSEIF(IOP2 .EQ. 2)THEN
  269.          DO 20 K=1,NMAI
  270.             IPT1=MCHELM.IMACHE(K)
  271.             SEGACT IPT1
  272.             NBNN=IPT1.NUM(/1)
  273.             NTOTV=NBNN
  274.             NBELEM=IPT1.NUM(/2)
  275.             MCHAM1=MCHELM.ICHAML(K)
  276.             SEGACT MCHAM1
  277.             MELVX=MCHAM1.IELVAL(1)
  278.             MELVY=MCHAM1.IELVAL(2)
  279.             SEGACT MELVX
  280.             SEGACT MELVY
  281.             IF (IDIM .EQ. 3) THEN
  282.                MELVZ=MCHAM1.IELVAL(3)
  283.                SEGACT MELVZ
  284.             ENDIF
  285.             SEGDES MCHAM1
  286.             DO NLELEM = 1, NBELEM
  287.                NGCE = IPT1.NUM(NBNN,NLELEM)
  288.                NLCE = MLECEN.LECT(NGCE)
  289.                DO NVOI = 1, NTOTV
  290.                   NGCV = IPT1.NUM(NVOI,NLELEM)
  291.                   NLCV = MLECEN.LECT(NGCV)
  292.                   IF(NLCV .EQ. 0)THEN
  293.                      NLCL = MLENCL.LECT(NGCV)
  294. C
  295. C************* NGCV n'est pas un point centre, mais un point face.
  296. C
  297.                      IF(NLCL .NE. 0)THEN
  298. C
  299. C************************* NGCV est un point du SPG de C.L.
  300. C
  301.                         DO I3 = 1, NCOMP
  302.                            VCOMP(I3) =  MPOVCL.VPOCHA(NLCL,I3)
  303.                         ENDDO
  304.                      ELSE
  305.                         NLCF = MLEFAC.LECT(NGCV)
  306.                         IF (IDIM .EQ.2) THEN
  307.                            VCOMP(1) = MPOCHP.VPOCHA(NLCE,1)
  308.                            VCOMP(2) = MPOCHP.VPOCHA(NLCE,2)
  309.                            NX = MPONOR.VPOCHA(NLCF,1)
  310.                            NY = MPONOR.VPOCHA(NLCF,2)
  311.                            VN = VCOMP(1)*NX+VCOMP(2)*NY
  312.                            VCOMP(1) = VCOMP(1)-2*NX*VN
  313.                            VCOMP(2) = VCOMP(2)-2*NY*VN
  314.                         ELSE
  315.                            VCOMP(1) = MPOCHP.VPOCHA(NLCE,1)
  316.                            VCOMP(2) = MPOCHP.VPOCHA(NLCE,2)
  317.                            VCOMP(3) = MPOCHP.VPOCHA(NLCE,3)
  318.                            NX = MPONOR.VPOCHA(NLCF,1)
  319.                            NY = MPONOR.VPOCHA(NLCF,2)
  320.                            NZ = MPONOR.VPOCHA(NLCF,3)
  321.                            VN = VCOMP(1)*NX+VCOMP(2)*NY+VCOMP(3)*NZ
  322.                            VCOMP(1) = VCOMP(1)-2*NX*VN
  323.                            VCOMP(2) = VCOMP(2)-2*NY*VN
  324.                            VCOMP(3) = VCOMP(3)-2*NZ*VN
  325.                         ENDIF
  326.                      ENDIF
  327.                      DO I3 = 1, NCOMP
  328.                         ICOM = IDIM*(I3 -1)+1
  329.                         VAL = VCOMP(I3)
  330.                         MPOGRA.VPOCHA(NLCE,ICOM)= MPOGRA.VPOCHA(NLCE
  331.      $                       ,ICOM)+VAL * MELVX.VELCHE(NVOI
  332.      $                       ,NLELEM)
  333.                         MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE
  334.      $                       ,ICOM+1) +VAL * MELVY.VELCHE(NVOI
  335.      $                       ,NLELEM)
  336.                         IF (IDIM .EQ. 3) THEN
  337.                            MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE
  338.      $                          ,ICOM+2) +VAL * MELVZ.VELCHE(NVOI
  339.      $                          ,NLELEM)
  340.                         ENDIF
  341.                         MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3
  342.      $                       ),VAL)
  343.                         MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3
  344.      $                       ),VAL)
  345.                      ENDDO
  346.                   ELSE
  347. C
  348. C********** NGCV est un point centre
  349. C
  350.                      DO I3 = 1, NCOMP
  351.                         ICOM = IDIM*(I3 -1)+1
  352.                         VAL = MPOCHP.VPOCHA(NLCV,I3)
  353.                         MPOGRA.VPOCHA(NLCE,ICOM)=MPOGRA.VPOCHA(NLCE,ICOM
  354.      $                       )+VAL * MELVX.VELCHE(NVOI,NLELEM)
  355.                         MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE
  356.      $                       ,ICOM+1)+VAL * MELVY.VELCHE(NVOI,NLELEM)
  357.                         IF (IDIM .EQ. 3) THEN
  358.                            MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE
  359.      $                          ,ICOM+2)+VAL * MELVZ.VELCHE(NVOI,NLELEM
  360.      $                          )
  361.                         ENDIF
  362.                         MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3
  363.      $                       ),VAL)
  364.                         MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3
  365.      $                       ),VAL)
  366.                      ENDDO
  367.                   ENDIF
  368.                ENDDO
  369. C
  370. C******* Fin boucle sur les voisins de NGCE
  371. C
  372.             ENDDO
  373.             SEGDES MELVX
  374.             SEGDES MELVY
  375.             IF (IDIM .EQ. 3) SEGDES MELVZ
  376.             SEGDES IPT1
  377.  20      CONTINUE
  378.       ENDIF
  379.       SEGDES MPONOR
  380.       SEGDES MCHELM
  381.       RETURN
  382.       END
  383.  
  384.  
  385.  
  386.  
  387.  
  388.  

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