Esope source files

   1 : C KSOF      SOURCE    CHAT      05/01/13    01:08:43     5004           
   2 :       SUBROUTINE KSOF
   3 : C*************************************************************************
   4 : C
   5 : C  Objet   : Change un champoint VECT SOMMET en SCAL FACE
   6 : C  Syntaxe : VFACE = KSOF VSOMMET TABDOM ;
   7 : C
   8 : C*************************************************************************
   9 :       IMPLICIT INTEGER(I-N)
  10 :       IMPLICIT REAL*8 (A-H,O-Z)
  11 : 
  12 : -INC CCOPTIO
  13 : -INC SMCOORD
  14 : -INC SMLENTI
  15 :       POINTEUR LEF.MLENTI
  16 : -INC SMTABLE
  17 :       POINTEUR MTABD.MTABLE
  18 : -INC SMCHPOI
  19 :       POINTEUR MCHPF.MCHPOI,MCHPN.MCHPOI,MPOVF.MPOVAL,MPOVN.MPOVAL
  20 : -INC SMELEME
  21 :       POINTEUR MELEMF.MELEME,MELEMP.MELEME
  22 : 
  23 :       PARAMETER (NTB=1)
  24 :       CHARACTER*8 LTAB(NTB),TYPE,TYPC
  25 :       DIMENSION KTAB(NTB)
  26 :       DATA LTAB/'DOMAINE '/
  27 : C***
  28 :       IAXI=0
  29 :       IF(IFOMOD.EQ.0)IAXI=2
  30 : 
  31 :       NTO=NTB
  31 :       CALL LITABS(LTAB,KTAB,NTB,NTO,IRET)
  32 :       IF(IRET.EQ.0)RETURN
  33 :       MTABD=KTAB(1)
  34 :       SEGACT MTABD
  35 : 
  35 :       CALL LIROBJ('CHPOINT ',MCHPOI,1,IRET)
  36 :       IF(IRET.EQ.0)THEN
  37 :       WRITE(6,*)' On attend un CHPOINT'
  38 :       RETURN
  39 :       ENDIF
  39 :       CALL LICHT(MCHPOI,MPOVAL,TYPC,IGEOS)
  40 :       IF(TYPC.NE.'SOMMET')THEN
  41 :       WRITE(6,*)'On attend un CHAMPOINT SOMMET'
  42 :       SEGDES MCHPOI,MPOVAL
  43 :       RETURN
  44 :       ENDIF
  45 :       NC=VPOCHA(/2)
  46 :       IF(NC.NE.IDIM)THEN
  47 :       WRITE(6,*)'On attend un CHAMPOINT VECT SOMMET'
  48 :       SEGDES MCHPOI,MPOVAL
  49 :       RETURN
  50 :       ENDIF
  50 :       CALL KRIPAD(IGEOS,MLENTI)
  51 : 
  51 :       CALL LEKTAB(MTABD,'FACE',MELEMF)
  52 :       IF(MELEMF.EQ.0)GO TO 90
  52 : C?    CALL KRIPAD(MELEMF,MLENT1)
  53 :       SEGACT MELEMF
  54 :       NBF=MELEMF.NUM(/2)
  54 :       CALL RSETXI(LECT,MELEMF.NUM,NBF)
  55 : 
  55 :       CALL LEKTAB(MTABD,'FACEP',MELEMP)
  56 :       IF(MELEMP.EQ.0)GO TO 90
  56 :       CALL LEKTAB(MTABD,'XXNORMAF',MCHPN)
  57 :       IF(MCHPN.EQ.0)GO TO 90
  57 : C     CALL KRIPAD(MELEMP,MLENTI)
  58 : 
  58 :       CALL LICHT(MCHPN,MPOVN,TYPE,IGEOM)
  59 :       TYPE='FACE'
  60 :       NC=1
  60 :       CALL CRCHPT(TYPE,MELEMF,NC,MCHPF)
  60 :       CALL LICHT(MCHPF,MPOVF,TYPE,IGEOM)
  61 : 
  62 : 
  63 :       SEGACT MELEMP
  64 :       NBSOUS=MELEMP.LISOUS(/1)
  65 :       IF(NBSOUS.EQ.0)NBSOUS=1
  66 :       NF=0
  67 :       DO 1 KS=1,NBSOUS
  68 :       IPT1=MELEMP
  69 :       IF(NBSOUS.NE.1)IPT1=MELEMP.LISOUS(KS)
  70 :       SEGACT IPT1
  71 :       NS=IPT1.NUM(/1)-1
  72 :       NEL=IPT1.NUM(/2)
  73 : 
  74 :       IF(IAXI.EQ.0)THEN
  75 :       DO 2 K=1,NEL
  76 : C?    NF=NF+1
  77 : C?    NF=MLENT1.LECT(IPT1.NUM(NS+1,K))
  78 :       NF=       LECT(IPT1.NUM(NS+1,K))
  79 : C     write(6,*)' 0 NF=',nf
  80 : C     NF=IPT1.NUM(2,K)
  81 : C     NF=IPADL(NF)
  82 :       VI=0.D0
  83 :       DO 22 IS=1,NS
  84 :       N1=IPT1.NUM(IS,K)
  85 :       N1=LECT(N1)
  86 :       DO 21 I=1,IDIM
  87 :       VI=VI+VPOCHA(N1,I)*MPOVN.VPOCHA(NF,I)
  88 :  21   CONTINUE
  89 :  22   CONTINUE
  90 : C     write(6,*)' NF=',nf
  91 : C     write(6,*)'Normale : ',MPOVN.VPOCHA(NF,1),MPOVN.VPOCHA(NF,2)
  92 :       MPOVF.VPOCHA(NF,1)=VI/FLOAT(NS)
  93 : C     write(6,*)'VI/F=',MPOVF.VPOCHA(NF,1)
  94 :  2    CONTINUE
  95 : 
  96 :       ELSEIF(IAXI.EQ.2)THEN
  97 :       f23=2.D0/3.D0
  98 :       DO 3 K=1,NEL
  99 : C?      NF=NF+1
 100 : C?    NF=MLENT1.LECT(IPT1.NUM(NS+1,K))
 101 :       NF=       LECT(IPT1.NUM(NS+1,K))
 102 : C     write(6,*)' 2 NF=',nf
 103 :       N1=IPT1.NUM(1,K)
 104 :       R1=XCOOR((N1-1)*3+1)
 105 :       N1=LECT(N1)
 106 :       N2=IPT1.NUM(2,K)
 107 :       R2=XCOOR((N2-1)*3+1)
 108 :       N2=LECT(N2)
 109 :       VX1=VPOCHA(N1,1)
 110 :       VX2=VPOCHA(N2,1)
 111 :       VY1=VPOCHA(N1,2)
 112 :       VY2=VPOCHA(N2,2)
 113 :       DN = ABS(MPOVN.VPOCHA(NF,2))
 114 :       IF(DN.GT.1.D-6)THEN
 115 :          VN1=VX1*MPOVN.VPOCHA(NF,1)+VY1*MPOVN.VPOCHA(NF,2)
 116 :          VN2=VX2*MPOVN.VPOCHA(NF,1)+VY2*MPOVN.VPOCHA(NF,2)
 117 :          DR=R2-R1
 118 :          VN=VN1*R2-VN2*R1+F23*(VN2-VN1)/DR*(R2**3-R1**3)/(R2+R1)
 119 :          MPOVF.VPOCHA(NF,1)=VN/DR
 120 :       ELSE
 121 : C calcul simplifier pour les faces casi verticale (sinon X/0.)
 122 :          VX=(VX1+VX2)*0.5D0
 123 :          VY=(VY1+VY2)*0.5D0
 124 :          MPOVF.VPOCHA(NF,1)=VX*MPOVN.VPOCHA(NF,1)+VY*MPOVN.VPOCHA(NF,2)
 125 :       ENDIF
 126 :  3    CONTINUE
 127 :       ENDIF
 128 : 
 129 :       SEGDES IPT1
 130 :  1    CONTINUE
 131 :       SEGDES MELEMP,MELEMF
 132 : 
 133 :       SEGSUP MLENTI
 133 :       CALL ECROBJ('CHPOINT ',MCHPF)
 134 :       SEGDES MCHPF,MPOVF,MTABD
 135 :       SEGDES MCHPOI,MPOVAL
 136 :       SEGDES MCHPN,MPOVN
 137 : 
 138 :       RETURN
 139 : 
 140 :  90   CONTINUE
 141 :       WRITE(6,*)' Retour anormal de KSOF'
 142 :       RETURN
 143 : 
 144 :       END
 145 :  
 146 :  
 147 :