Sources Esope | Cast3M

ksof
C KSOF      SOURCE    CB215821  20/11/25    13:33:09     10792          
      SUBROUTINE KSOF
C*************************************************************************
C
C  Objet   : Change un champoint VECT SOMMET en SCAL FACE
C  Syntaxe : VFACE = KSOF VSOMMET TABDOM ;
C
C*************************************************************************
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
-INC SMLENTI
      POINTEUR LEF.MLENTI
-INC SMTABLE
      POINTEUR MTABD.MTABLE
-INC SMCHPOI
      POINTEUR MCHPF.MCHPOI,MCHPN.MCHPOI,MPOVF.MPOVAL,MPOVN.MPOVAL
-INC SMELEME
      POINTEUR MELEMF.MELEME,MELEMP.MELEME
 
      PARAMETER (NTB=1)
      CHARACTER*8 LTAB(NTB),TYPE,TYPC
      DIMENSION KTAB(NTB)
      DATA LTAB/'DOMAINE '/
C***
      IAXI=0
      IF(IFOMOD.EQ.0)IAXI=2
 
      NTO=NTB
      CALL LITABS(LTAB,KTAB,NTB,NTO,IRET)
      IF(IRET.EQ.0)RETURN
      MTABD=KTAB(1)
      SEGACT MTABD
 
      CALL LIROBJ('CHPOINT ',MCHPOI,1,IRET)
      IF(IRET.EQ.0)THEN
      WRITE(6,*)' On attend un CHPOINT'
      RETURN
      ENDIF
      CALL LICHT(MCHPOI,MPOVAL,TYPC,IGEOS)
      IF(TYPC.NE.'SOMMET')THEN
      WRITE(6,*)'On attend un CHAMPOINT SOMMET'
      SEGDES MCHPOI,MPOVAL
      RETURN
      ENDIF
      NC=VPOCHA(/2)
      IF(NC.NE.IDIM)THEN
      WRITE(6,*)'On attend un CHAMPOINT VECT SOMMET'
      SEGDES MCHPOI,MPOVAL
      RETURN
      ENDIF
      CALL KRIPAD(IGEOS,MLENTI)
 
      CALL LEKTAB(MTABD,'FACE',MELEMF)
      IF(MELEMF.EQ.0)GO TO 90
C?    CALL KRIPAD(MELEMF,MLENT1)
      SEGACT MELEMF
      NBF=MELEMF.NUM(/2)
      CALL RSETXI(LECT,MELEMF.NUM,NBF)
 
      CALL LEKTAB(MTABD,'FACEP',MELEMP)
      IF(MELEMP.EQ.0)GO TO 90
      CALL LEKTAB(MTABD,'XXNORMAF',MCHPN)
      IF(MCHPN.EQ.0)GO TO 90
C     CALL KRIPAD(MELEMP,MLENTI)
 
      CALL LICHT(MCHPN,MPOVN,TYPE,IGEOM)
      TYPE='FACE'
      NC=1
      CALL CRCHPT(TYPE,MELEMF,NC,MCHPF)
      CALL LICHT(MCHPF,MPOVF,TYPE,IGEOM)
 
 
      SEGACT MELEMP
      NBSOUS=MELEMP.LISOUS(/1)
      IF(NBSOUS.EQ.0)NBSOUS=1
      NF=0
      DO 1 KS=1,NBSOUS
      IPT1=MELEMP
      IF(NBSOUS.NE.1)IPT1=MELEMP.LISOUS(KS)
      SEGACT IPT1
      NS=IPT1.NUM(/1)-1
      NEL=IPT1.NUM(/2)
 
      IF(IAXI.EQ.0)THEN
      DO 2 K=1,NEL
C?    NF=NF+1
C?    NF=MLENT1.LECT(IPT1.NUM(NS+1,K))
      NF=       LECT(IPT1.NUM(NS+1,K))
C     write(6,*)' 0 NF=',nf
C     NF=IPT1.NUM(2,K)
C     NF=IPADL(NF)
      VI=0.D0
      DO 22 IS=1,NS
      N1=IPT1.NUM(IS,K)
      N1=LECT(N1)
      DO 21 I=1,IDIM
      VI=VI+VPOCHA(N1,I)*MPOVN.VPOCHA(NF,I)
 21   CONTINUE
 22   CONTINUE
C     write(6,*)' NF=',nf
C     write(6,*)'Normale : ',MPOVN.VPOCHA(NF,1),MPOVN.VPOCHA(NF,2)
      MPOVF.VPOCHA(NF,1)=VI/FLOAT(NS)
C     write(6,*)'VI/F=',MPOVF.VPOCHA(NF,1)
 2    CONTINUE
 
      ELSEIF(IAXI.EQ.2)THEN
      f23=2.D0/3.D0
      DO 3 K=1,NEL
C?      NF=NF+1
C?    NF=MLENT1.LECT(IPT1.NUM(NS+1,K))
      NF=       LECT(IPT1.NUM(NS+1,K))
C     write(6,*)' 2 NF=',nf
      N1=IPT1.NUM(1,K)
      R1=XCOOR((N1-1)*3+1)
      N1=LECT(N1)
      N2=IPT1.NUM(2,K)
      R2=XCOOR((N2-1)*3+1)
      N2=LECT(N2)
      VX1=VPOCHA(N1,1)
      VX2=VPOCHA(N2,1)
      VY1=VPOCHA(N1,2)
      VY2=VPOCHA(N2,2)
      DN = ABS(MPOVN.VPOCHA(NF,2))
      IF(DN.GT.1.D-6)THEN
         VN1=VX1*MPOVN.VPOCHA(NF,1)+VY1*MPOVN.VPOCHA(NF,2)
         VN2=VX2*MPOVN.VPOCHA(NF,1)+VY2*MPOVN.VPOCHA(NF,2)
         DR=R2-R1
         VN=VN1*R2-VN2*R1+F23*(VN2-VN1)/DR*(R2**3-R1**3)/(R2+R1)
         MPOVF.VPOCHA(NF,1)=VN/DR
      ELSE
C calcul simplifier pour les faces casi verticale (sinon X/0.)
         VX=(VX1+VX2)*0.5D0
         VY=(VY1+VY2)*0.5D0
         MPOVF.VPOCHA(NF,1)=VX*MPOVN.VPOCHA(NF,1)+VY*MPOVN.VPOCHA(NF,2)
      ENDIF
 3    CONTINUE
      ENDIF
 
      SEGDES IPT1
 1    CONTINUE
      SEGDES MELEMP,MELEMF
 
      SEGSUP MLENTI
      CALL ECROBJ('CHPOINT ',MCHPF)
      SEGDES MCHPF,MPOVF,MTABD
      SEGDES MCHPOI,MPOVAL
      SEGDES MCHPN,MPOVN
 
      RETURN
 
 90   CONTINUE
      WRITE(6,*)' Retour anormal de KSOF'
      RETURN
 
      END