C PENTE2 SOURCE CB215821 20/11/25 13:35:38 10792 SUBROUTINE PENTE2(IOP2,INORM,MLECEN,MLEFAC, & MLENCL,IMCHAM, & NCOMP,MPOCHP,MPOVCL,MPOGRA, & MPOMIN,MPOMAX) C C************************************************************************ C C PROJET : CASTEM 2000 C C NOM : PENTE2 C C DESCRIPTION : Cette subroutine est appellée par la subroutine C PENTE1 (calcul du gradient d'un CHPOINT de type C CENTRE) C Elle contient la partie du calcul du gradient. C C LANGUAGE : ESOPE 2000 avec extensions CISI C C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF C AUTEUR (modif.) : R. MOREL, DRN/DMT/SEMT/TTMF C C************************************************************************ C C C APPELES (Outils) : LICHT, ERREUR C C C************************************************************************ C C ENTREES : IOP2 : Entier indiquant le type de reconstruction (cf PENT) C C INORM : Pointeur sur le CHPOINT des normales aux faces C C MLECEN : table numerotation global/local pour les C points CENTRE (segment MLENTI) C C MLEFAC : table numerotation global/local pour les C points FACE C C MLENCL : table numerotation global/local pour les C SPG du CHPOINT de C.L. ou C table de 0 si ce CHPOINT n'existe pas C C IMCHAM : MCHAML de type 'GRADGEO' C C NCOMP : Nombre de composantes du CHPOINT dont C on veut le gradient C C MPOCHP : MPOVAL du CHPOINT dont on veut le gradient C C MPOVCL : MPOVAL du CHPOINT de C.L. C C SORTIES : MPOGRA : MPOVAL du gradient C C MPOMAX : MPOVAL du max sur le stencil C C MPOMIN : MPOVAL du min sur le stencil C C C************************************************************************ C C HISTORIQUE (Anomalies et modifications éventuelles) C C HISTORIQUE : Cree le 4-6-1998 C C HISTORIQUE : Modifie le 20-10-1998 pour extension 3D C Modifie le 28-04-2000 pour reconstruction quadatique C exacte (A. Beccantini) C HISTORIQUE : Modifie le 05.07.2003: elimination reconstruction C quadratique exacte. C Algorithme de BArth-Jespersen remplacé par une C reconstruction lineaire exacte qui se base sur C LSM C************************************************************************ C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC SMCHAML -INC SMELEME -INC SMCHPOI -INC SMLENTI C POINTEUR INORM.MCHPOI POINTEUR MLENCL.MLENTI, MLECEN.MLENTI, MLEFAC.MLENTI POINTEUR MPOMAX.MPOVAL, MPOMIN.MPOVAL, MPONOR.MPOVAL, & MPOVCL.MPOVAL, MPOCHP.MPOVAL, MPOGRA.MPOVAL POINTEUR MELVX.MELVAL, MELVY.MELVAL, MELVZ.MELVAL, & MELVXX.MELVAL, MELVYY.MELVAL, MELVZZ.MELVAL, & MELVXY.MELVAL, MELVXZ.MELVAL, MELVYZ.MELVAL C INTEGER IOP2,NCOMP,NLELEM,NVOI,I3,IMCHAM & ,NGCE,NLCE,NTOTV,NGCV,NLCV,ICOM,IGEOM,K,NBELEM,NBNN & ,NLCF,NLCL,NMAI REAL*8 VAL,VCOMP(3),VN REAL*8 NX,NY,NZ CHARACTER*8 TYPE C C**** Segments déjà activés C C MPOVCL C MLENCL C MLECEN C MLEFAC C MPOCHP C MPOGRA*MOD C MPOMIN*MOD C MPOMAX*MOD C C**** N.B. MPOMAX.VPOCHA, MPOMIN.VPOCHA déjà initialisé: C C MPOMAX.VPOCHA = MPOCHP.VPOCHA C MPOMIN.VPOCHA = MPOCHP.VPOCHA C CALL LICHT(INORM,MPONOR,TYPE,IGEOM) C C**** En LICHT C SEGACT*MOD MPONOR C MPONOR: pointeur sur les normales aux faces C MCHELM = IMCHAM SEGACT MCHELM NMAI=MCHELM.ICHAML(/1) C C********************************************* C Premier cas : on traite un champ scalaire en C reconstruction non-quadratique C********************************************* C C IF(IOP2.EQ.1) THEN C C**** N.B: IOP2 = 1 -> EULESCAL C DO 10 K=1,NMAI C C********** Boucle sur les differents maillages elementaires C IPT1 = maillage duale C IPT1=MCHELM.IMACHE(K) SEGACT IPT1 NBNN=IPT1.NUM(/1) NTOTV=NBNN C C********** NTOTV = Nombre de voisins C NBELEM = IPT1.NUM(/2) MCHAM1 = MCHELM.ICHAML(K) SEGACT MCHAM1 MELVX = MCHAM1.IELVAL(1) MELVY = MCHAM1.IELVAL(2) SEGACT MELVX SEGACT MELVY IF (IDIM .EQ. 3) THEN MELVZ = MCHAM1.IELVAL(3) SEGACT MELVZ ENDIF SEGDES MCHAM1 DO NLELEM = 1, NBELEM, 1 C C************* Boucle sur les ELTs C NLELEM = Numero local d'ELT dans le maillage 1 C C NGCE = numero global centre d'elt C NLCE = numero local centre d'elt C C MLECEN = table numerotation global/local pour les C points CENTRE (segment MLENTI) C C MLEFAC = table numerotation global/local pour les C points FACE C C MLENCL = table numerotation global/local pour les C SPG du CHPOINT de C.L. ou C table de 0 si ce CHPOINT n'existe pas C NGCE = IPT1.NUM(NBNN,NLELEM) NLCE = MLECEN.LECT(NGCE) DO NVOI = 1, NTOTV C C**************** Boucle sur les voisins de NGCE C C NGCV = Numero global voisin C NLCV = " local de NGCV dans le maillage CENTRE C NLCL = " local de NGCV dans le SPG du CHPOINT C des C.L. C NGCV = IPT1.NUM(NVOI,NLELEM) NLCV = MLECEN.LECT(NGCV) IF(NLCV.EQ.0)THEN NLCL = MLENCL.LECT(NGCV) C C**************** NGCV n'est pas un point centre, mais un point face, C et on est sur les bords C DO I3 = 1, NCOMP C C********************** I3 : numero de composante du CHPOI C ICOM : numero de composante du gradient C ICOM = IDIM*(I3 -1)+1 IF (NLCL .EQ. 0) THEN C C************************* NGCV n'est pas un point du SPG de C.L. C VAL = MPOCHP.VPOCHA(NLCE,I3) ELSE C C************************* NGCV est un point du SPG de C.L. C VAL = MPOVCL.VPOCHA(NLCL,I3) ENDIF MPOGRA.VPOCHA(NLCE,ICOM)= MPOGRA.VPOCHA(NLCE $ ,ICOM)+VAL * MELVX.VELCHE(NVOI,NLELEM) MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE $ ,ICOM+1) +VAL * MELVY.VELCHE(NVOI,NLELEM) IF (IDIM .EQ. 3) THEN MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE $ ,ICOM+2) +VAL * MELVZ.VELCHE(NVOI $ ,NLELEM) ENDIF MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3 $ ),VAL) MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3 $ ),VAL) ENDDO ELSE C C********** NGCV est un point centre C DO I3 = 1, NCOMP ICOM = IDIM*(I3 -1)+1 VAL = MPOCHP.VPOCHA(NLCV,I3) MPOGRA.VPOCHA(NLCE,ICOM)=MPOGRA.VPOCHA(NLCE,ICOM $ )+VAL * MELVX.VELCHE(NVOI,NLELEM) MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE $ ,ICOM+1)+VAL * MELVY.VELCHE(NVOI,NLELEM) IF (IDIM .EQ. 3) THEN MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE $ ,ICOM+2)+VAL * MELVZ.VELCHE(NVOI,NLELEM $ ) ENDIF MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3 $ ),VAL) MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3 $ ),VAL) ENDDO ENDIF ENDDO C C******* Fin boucle sur les voisins de NGCE C ENDDO SEGDES MELVX SEGDES MELVY IF (IDIM .EQ. 3) SEGDES MELVZ SEGDES IPT1 10 CONTINUE C C********************************************************** C Deuxieme cas : on traite un champ vectoriel (IOP2=2) C********************************************************** C C ELSEIF(IOP2 .EQ. 2)THEN DO 20 K=1,NMAI IPT1=MCHELM.IMACHE(K) SEGACT IPT1 NBNN=IPT1.NUM(/1) NTOTV=NBNN NBELEM=IPT1.NUM(/2) MCHAM1=MCHELM.ICHAML(K) SEGACT MCHAM1 MELVX=MCHAM1.IELVAL(1) MELVY=MCHAM1.IELVAL(2) SEGACT MELVX SEGACT MELVY IF (IDIM .EQ. 3) THEN MELVZ=MCHAM1.IELVAL(3) SEGACT MELVZ ENDIF SEGDES MCHAM1 DO NLELEM = 1, NBELEM NGCE = IPT1.NUM(NBNN,NLELEM) NLCE = MLECEN.LECT(NGCE) DO NVOI = 1, NTOTV NGCV = IPT1.NUM(NVOI,NLELEM) NLCV = MLECEN.LECT(NGCV) IF(NLCV .EQ. 0)THEN NLCL = MLENCL.LECT(NGCV) C C************* NGCV n'est pas un point centre, mais un point face. C IF(NLCL .NE. 0)THEN C C************************* NGCV est un point du SPG de C.L. C DO I3 = 1, NCOMP VCOMP(I3) = MPOVCL.VPOCHA(NLCL,I3) ENDDO ELSE NLCF = MLEFAC.LECT(NGCV) IF (IDIM .EQ.2) THEN VCOMP(1) = MPOCHP.VPOCHA(NLCE,1) VCOMP(2) = MPOCHP.VPOCHA(NLCE,2) NX = MPONOR.VPOCHA(NLCF,1) NY = MPONOR.VPOCHA(NLCF,2) VN = VCOMP(1)*NX+VCOMP(2)*NY VCOMP(1) = VCOMP(1)-2*NX*VN VCOMP(2) = VCOMP(2)-2*NY*VN ELSE VCOMP(1) = MPOCHP.VPOCHA(NLCE,1) VCOMP(2) = MPOCHP.VPOCHA(NLCE,2) VCOMP(3) = MPOCHP.VPOCHA(NLCE,3) NX = MPONOR.VPOCHA(NLCF,1) NY = MPONOR.VPOCHA(NLCF,2) NZ = MPONOR.VPOCHA(NLCF,3) VN = VCOMP(1)*NX+VCOMP(2)*NY+VCOMP(3)*NZ VCOMP(1) = VCOMP(1)-2*NX*VN VCOMP(2) = VCOMP(2)-2*NY*VN VCOMP(3) = VCOMP(3)-2*NZ*VN ENDIF ENDIF DO I3 = 1, NCOMP ICOM = IDIM*(I3 -1)+1 VAL = VCOMP(I3) MPOGRA.VPOCHA(NLCE,ICOM)= MPOGRA.VPOCHA(NLCE $ ,ICOM)+VAL * MELVX.VELCHE(NVOI $ ,NLELEM) MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE $ ,ICOM+1) +VAL * MELVY.VELCHE(NVOI $ ,NLELEM) IF (IDIM .EQ. 3) THEN MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE $ ,ICOM+2) +VAL * MELVZ.VELCHE(NVOI $ ,NLELEM) ENDIF MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3 $ ),VAL) MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3 $ ),VAL) ENDDO ELSE C C********** NGCV est un point centre C DO I3 = 1, NCOMP ICOM = IDIM*(I3 -1)+1 VAL = MPOCHP.VPOCHA(NLCV,I3) MPOGRA.VPOCHA(NLCE,ICOM)=MPOGRA.VPOCHA(NLCE,ICOM $ )+VAL * MELVX.VELCHE(NVOI,NLELEM) MPOGRA.VPOCHA(NLCE,ICOM+1)=MPOGRA.VPOCHA(NLCE $ ,ICOM+1)+VAL * MELVY.VELCHE(NVOI,NLELEM) IF (IDIM .EQ. 3) THEN MPOGRA.VPOCHA(NLCE,ICOM+2)=MPOGRA.VPOCHA(NLCE $ ,ICOM+2)+VAL * MELVZ.VELCHE(NVOI,NLELEM $ ) ENDIF MPOMAX.VPOCHA(NLCE,I3)=MAX(MPOMAX.VPOCHA(NLCE,I3 $ ),VAL) MPOMIN.VPOCHA(NLCE,I3)=MIN(MPOMIN.VPOCHA(NLCE,I3 $ ),VAL) ENDDO ENDIF ENDDO C C******* Fin boucle sur les voisins de NGCE C ENDDO SEGDES MELVX SEGDES MELVY IF (IDIM .EQ. 3) SEGDES MELVZ SEGDES IPT1 20 CONTINUE ENDIF SEGDES MPONOR SEGDES MCHELM RETURN END