Télécharger fptau.eso

Retour à la liste

fptau
C FPTAU     SOURCE    PV        08/09/11    21:15:54     6150
      SUBROUTINE FPTAU(IDIM,NP,NBEL,NUM,IPADT,XCOOR,
     & U,NPTI,NC,IPADS,NPTS,ANUT,AMUEF,
     & YP,AN,RO,IKR,AMU,IKM,UET,YPS,ITYPEL,TKE,NBINC)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C.......................................................................
C  Ut(NP,NBEL)       module vitesse tangente = |U - Un|
C  Un(NP,NBEL)       module vitesse normale  = Un
C  YP(NP)            distance la paroi
C  UET(*)            uet aux noeuds frontieres
C  AN(IDIM,NP)  normale
C  ANN               module de la normale
C  NUM(NP)           connectivite element appuyé
C  NMS(NS)           connectivite relative (frontiere)
C  NMY(NY)           connectivite complementaire
C  RO(*)             densite
C  AMUEF(NP,NBEL)    Viscosite effective aux noeuds/elt
C
C
C1/ ce qu'il faudrait faire c'est attribuer une normale pour chaque point
C   de l'element appuye -> c'est fait
C.......................................................................
 
-INC CCREEL
      DIMENSION LKKE(1000)
      real*8 ann(3),Xa(9),Ya(9),Za(9)
      DIMENSION XCOOR(*),XYZ(3,27)
      DIMENSION NUM(NP,NBEL),NMS(9),NMY(27)
      DIMENSION U(NPTI,NC),IPADT(*),VS(27,3),IPADS(*)
      DIMENSION ANUT(NPTI),AMUEF(NP,NBEL)
      DIMENSION YP(NP,NBEL),AN(IDIM,NP,NBEL)
      DIMENSION UT(27),UN(27)
      DIMENSION RO(*),AMU(*)
      DIMENSION UET(NPTS),YPS(NPTS),TKE(NPTS,*)
      DIMENSION NUMFA(27)
      DATA EPS/1.D-4/,NBIT/100/
c-----------------------------------------------------------------------
c------- Caracterisation des elements ----------------------------------
c-----------------------------------------------------------------------
      DIMENSION NMF9(3,4),N1F9(3,4),N2F9(3,4)
      DATA NMF9/1,2,3,3,4,5,5,6,7,7,8,1/
      DATA N1F9/8,9,4,2,9,6,4,9,8,6,9,2/,N2F9/7,6,5,1,8,7,3,2,1,5,4,3/
c-----------------------------------------------------------------------
      DIMENSION NMF7(3,3),N1F7(3,3),N2F7(1,3)
      DATA NMF7/1,2,3,3,4,5,5,6,1/
      DATA N1F7/6,7,4,2,7,6,4,7,2/,N2F7/5,1,3/
c-----------------------------------------------------------------------
      DIMENSION NMF8(4,6),N1F8(4,6),N2F8(4,6)
      DATA NMF8/1,2,3,4,2,3,7,6,5,6,7,8,1,4,8,5,1,2,6,5,4,3,7,8/
c?    DATA N1F8/6,7,4,2,7,6,4,7,2/,N2F8/5,1,3/
c-----------------------------------------------------------------------
C CONSTANTES PHYSIQUES
      DATA CNU,AKER /0.09D0,0.41D0/
c-----------------------------------------------------------------------
c-----------------------------------------------------------------------
 
      SQCNU=SQRT(CNU)
 
c     write(6,*)' DEBUT FPTAU : NBINC=',NBINC
c     write(6,*)'NES=',nes,' NPG=',npg,' idim=',idim,' np=',np
c     write(6,*)'NPTS=',NPTS
 
      DEUPI=1.D0
cc    IF(IAXI.NE.0)DEUPI=2.D0*XPI
      IERC=0
      ERRMAX=0.D0
 
      KT=0
      DO 2 K=1,NBEL
c     write(6,*)' DEBUT BCL sur K ',K
      NS=0
      NY=0
      CALL INITI(NMS,NP,0)
c     CALL INITI(MMS,NP,0)
      CALL INITI(NMY,NP,0)
c     write(6,*)' ELEMENT K=',K
 
      DO 20 I=1,NP
      J=NUM(I,K)
      JI=IPADT(NUM(I,K))
      JM=(1-IKM)*(JI-1)+1
c     JS=IPADS(NUM(I,K))
c      write(6,*)'I=',I,' num(i=',NUM(I,K),' Ipads()=',
c    & IPADS(J)
       IF(IPADS(J).NE.0)THEN
       NS=NS+1
       NMS(NS)=I
c      MMS(NS)=I
c      write(6,*)'I=',I,' num(i=',NUM(I,K),' Ipads()=',
c    & IPADS(J),' NS=',NS,'  ',NMS(NS)
       ELSE
       NY=NY+1
       NMY(NY)=I
       ENDIF
        DO 10 N=1,IDIM
          XYZ(N,I) = XCOOR((J-1)*(IDIM+1)+N)
          VS(I,N)=U(JI,N)
  10    CONTINUE
c              UETK(I,K)=VET(JS)
c              UETK(I,K)=0.D0
c              ROS(I) = RO(JR)
c              AMUS(I)= AMU(JM)
  20     CONTINUE
c      write(6,*)'K=',K,' MMS=',(mms(ijk),ijk=1,nS)
 
c      write(6,*)'K=',K,' XYZ=',(XYZ(1,I),I=1,np)
c      write(6,*)'K=',K,' XYZ=',(XYZ(2,I),I=1,np)
c      write(6,*)'K=',K,' NP=',NP,' NS=',NS,' NY=',NY
c      write(6,*)'K=',K,' NUM ',(num(i,k),i=1,np)
c      write(6,*)'K=',K,' NMSO',(nms(i),i=1,nS)
c      write(6,*)'K=',K,' NMS ',(num(nms(i),k),i=1,nS)
c      write(6,*)'K=',K,' NMY ',(num(NMY(I),k),I=1,NY)
 
       IF(NY.EQ.1)GO TO 2
       IF(NS.EQ.1)GO TO 2
       IF(NS.LE.2.AND.IDIM.EQ.3)GO TO 2
       KT=KT+1
c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
c._._ On Calcule la normale au SEG2 (NMS) et YP
c._._ (La normale est la meme pour les deux pts)
      IF(IDIM.EQ.2.AND.NS.EQ.2)THEN
 
      Xa(1)=XYZ(1,NMS(1))
      Ya(1)=XYZ(2,NMS(1))
      Xa(2)=XYZ(1,NMS(2))
      Ya(2)=XYZ(2,NMS(2))
c     a=(Ya(1)-Ya(2))/(Xa(1)-Xa(2))
c     a=1.
      aa=Ya(1)-Ya(2)
      bb=Xa(1)-Xa(2)
      AN(1,NMS(1),K)=aa
      AN(2,NMS(1),K)=(-1.)*bb
      ann(1)=((aa*aa)+(bb*bb))**0.5D0
      AN(1,NMS(1),K)=AN(1,NMS(1),K)/ann(1)
      AN(2,NMS(1),K)=AN(2,NMS(1),K)/ann(1)
      DO 241 I=1,NY
      Xq=XYZ(1,NMY(I))
      Yq=XYZ(2,NMY(I))
      YP(NMY(I),K)=abs((aa*Xq)-(bb*Yq)+(bb*Ya(1))-(aa*Xa(1)))/ann(1)
 241  CONTINUE
c     la normale est la meme pour tous les pts (SEG2)
      DO 251 I=2,NS
       AN(1,NMS(I),K)=AN(1,NMS(1),K)
       AN(2,NMS(I),K)=AN(2,NMS(1),K)
 251  CONTINUE
      DO 252 I=1,NY
       AN(1,NMY(I),K)=AN(1,NMS(1),K)
       AN(2,NMY(I),K)=AN(2,NMS(1),K)
 252  CONTINUE
c._._ FIN On Calcule la normale au SEG2 (NMS) et YP
c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
 
c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
c._._ On Calcule la normale au SEG3 (NMS) et YP
c._._ (La normale n'est pas la meme pour tous les pts)
      ELSEIF(IDIM.EQ.2.AND.NS.EQ.3)THEN
 
C._._ Mais quelle est donc la face?
      call initi(numfa,np,0)
      do 2300 I=1,NS
      numfa(nms(I))=1
 2300 continue
      nfa=numfa(2)+numfa(4)*2+numfa(6)*3+numfa(8)*4
c     write(6,*)'NFA=',NFA
 
      Xa(1)=XYZ(1,NMF9(1,NFA))
      Ya(1)=XYZ(2,NMF9(1,NFA))
      Xa(2)=XYZ(1,NMF9(2,NFA))
      Ya(2)=XYZ(2,NMF9(2,NFA))
      Xa(3)=XYZ(1,NMF9(3,NFA))
      Ya(3)=XYZ(2,NMF9(3,NFA))
 
      a=(Ya(1)-Ya(2))/(Xa(1)-Xa(2))
      AN(1,NUM(NMF9(1,NFA),K),K)=a
      AN(2,NUM(NMF9(1,NFA),K),K)=-1.D0
      ann(1)=(a*a+1.D0)**0.5D0
      AN(1,NUM(NMF9(1,NFA),K),K)=AN(1,NUM(NMF9(1,NFA),K),K)/ann(1)
      AN(2,NUM(NMF9(1,NFA),K),K)=AN(2,NUM(NMF9(1,NFA),K),K)/ann(1)
 
      b=(Ya(1)-Ya(3))/(Xa(1)-Xa(3))
      AN(1,NUM(NMF9(2,NFA),K),K)=b
      AN(2,NUM(NMF9(2,NFA),K),K)=-1.D0
      ann(2)=(b*b+1.D0)**0.5D0
      AN(1,NUM(NMF9(2,NFA),K),K)=AN(1,NUM(NMF9(2,NFA),K),K)/ann(2)
      AN(2,NUM(NMF9(2,NFA),K),K)=AN(2,NUM(NMF9(2,NFA),K),K)/ann(2)
 
      c=(Ya(3)-Ya(2))/(Xa(3)-Xa(2))
      AN(1,NUM(NMF9(3,NFA),K),K)=c
      AN(2,NUM(NMF9(3,NFA),K),K)=-1.D0
      ann(3)=(c*c+1.D0)**0.5D0
      AN(1,NUM(NMF9(3,NFA),K),K)=AN(1,NUM(NMF9(3,NFA),K),K)/ann(3)
      AN(2,NUM(NMF9(3,NFA),K),K)=AN(2,NUM(NMF9(3,NFA),K),K)/ann(3)
 
c     write(6,*)'ann=',(ann(kk),kk=1,ns)
      DO 231 I=1,NS
      Xq=XYZ(1,N1F9(I,NFA))
      Yq=XYZ(2,N1F9(I,NFA))
      YP(N1F9(I,NFA),K)=abs(a*Xq-Yq+Ya(I)-a*Xa(I))/ann(I)
      Xq=XYZ(1,N2F9(I,NFA))
      Yq=XYZ(2,N2F9(I,NFA))
      YP(N2F9(I,NFA),K)=abs(a*Xq-Yq+Ya(I)-a*Xa(I))/ann(I)
 231  CONTINUE
 
c     write(6,*)'NMY=',(NMY(I),I=1,NY)
c     write(6,*)'YP=',(YP(NMY(I),K),I=1,NY)
c     write(6,*)'YP=',(YP(NMS(I),K),I=1,NS)
c     la normale n'est pas la meme pour tous les pts (SEG3)
      DO 232 I=1,NS
       AN(1,N1F9(I,NFA),K)=AN(1,NMS(I),K)
       AN(2,N1F9(I,NFA),K)=AN(2,NMS(I),K)
       AN(1,N2F9(I,NFA),K)=AN(1,NMS(I),K)
       AN(2,N2F9(I,NFA),K)=AN(2,NMS(I),K)
 232  CONTINUE
c._._ On Calcule la normale au SEG3 (NMS) et YP
c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
 
c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
c._._ On Calcule les normales aux CUB8,PRI6,TET4 et PYR5
c._._ (La normale est la meme pour tous les pts)
      ELSEIF(    ITYPEL.EQ.14
     &       .OR.ITYPEL.EQ.16
     &       .OR.ITYPEL.EQ.23
     &       .OR.ITYPEL.EQ.25 )THEN
 
      Xa(1)=XYZ(1,NMS(1))
      Ya(1)=XYZ(2,NMS(1))
      Za(1)=XYZ(3,NMS(1))
 
      Xa(2)=XYZ(1,NMS(2))
      Ya(2)=XYZ(2,NMS(2))
      Za(2)=XYZ(3,NMS(2))
 
      Xa(3)=XYZ(1,NMS(3))
      Ya(3)=XYZ(2,NMS(3))
      Za(3)=XYZ(3,NMS(3))
 
      AN(1,NMS(1),K)=
     & (Ya(2)-Ya(1))*(Za(3)-Za(1))-(Za(2)-Za(1))*(Ya(3)-Ya(1))
      AN(2,NMS(1),K)=
     & (Za(2)-Za(1))*(Xa(3)-Xa(1))-(Xa(2)-Xa(1))*(Za(3)-Za(1))
      AN(3,NMS(1),K)=
     & (Xa(2)-Xa(1))*(Za(3)-Za(1))-(Ya(2)-Ya(1))*(Xa(3)-Xa(1))
      ann(1)=(AN(1,NMS(1),K)*AN(1,NMS(1),K)
     &       +AN(2,NMS(1),K)*AN(2,NMS(1),K)
     &       +AN(3,NMS(1),K)*AN(3,NMS(1),K))**0.5D0
      AN(1,NMS(1),K)=AN(1,NMS(1),K)/ann(1)
      AN(2,NMS(1),K)=AN(2,NMS(1),K)/ann(1)
      AN(3,NMS(1),K)=AN(3,NMS(1),K)/ann(1)
 
      DO 431 I=2,NS
      AN(1,NMS(I),K)=AN(1,NMS(1),K)
      AN(2,NMS(I),K)=AN(2,NMS(1),K)
      AN(3,NMS(I),K)=AN(3,NMS(1),K)
 431  CONTINUE
 
c     write(6,*)'ann=',ann(1)
      DO 432 I=1,NY
      Xq=XYZ(1,NMY(I))
      Yq=XYZ(2,NMY(I))
      Zq=XYZ(3,NMY(I))
      YP(NMY(I),K)=ABS(
     &     (Xq-Xa(1))*AN(1,NMS(1),K)+(Yq-Ya(1))*AN(2,NMS(1),K)
     &                              +(Zq-Za(1))*AN(3,NMS(1),K))
c     write(6,*)'(Xq-Xa(1))=',(Xq-Xa(1)),' (Yq-Ya(1))=',(Yq-Ya(1)),
c    &' (Zq-Za(1))=',(Zq-Za(1))
 432  CONTINUE
 
c     write(6,*)'AN ',AN(1,NMS(1),K),AN(2,NMS(1),K),AN(3,NMS(1),K)
 
c     write(6,*)'NMY=',(NMY(I),I=1,NY)
c     write(6,*)'YP=',(YP(NMY(I),K),I=1,NY)
c     write(6,*)'YP=',(YP(NMS(I),K),I=1,NS)
c     la normale n'est pas la meme pour tous les pts (SEG3)
      DO 433 I=1,NY
       AN(1,NMY(I),K)=AN(1,NMS(1),K)
       AN(2,NMY(I),K)=AN(2,NMS(1),K)
       AN(3,NMY(I),K)=AN(3,NMS(1),K)
 433  CONTINUE
c._._ FIN Calcul des normales aux CUB8,PRI6,TET4 et PYR5
c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
 
c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
c._._ On Calcule les normales aux CU27,PR21,TE15 et PY19
c._._ (La normale est la meme pour tous les pts)
      ELSEIF(    ITYPEL.EQ.33
     &       .OR.ITYPEL.EQ.34
     &       .OR.ITYPEL.EQ.35
     &       .OR.ITYPEL.EQ.36 )THEN
 
      Xa(1)=XYZ(1,NMS(1))
      Ya(1)=XYZ(2,NMS(1))
      Za(1)=XYZ(3,NMS(1))
 
      Xa(2)=XYZ(1,NMS(2))
      Ya(2)=XYZ(2,NMS(2))
      Za(2)=XYZ(3,NMS(2))
 
      Xa(3)=XYZ(1,NMS(3))
      Ya(3)=XYZ(2,NMS(3))
      Za(3)=XYZ(3,NMS(3))
 
      Xa(4)=XYZ(1,NMS(NS))
      Ya(4)=XYZ(2,NMS(NS))
      Za(4)=XYZ(3,NMS(NS))
 
      ANx1=(Ya(2)-Ya(1))*(Za(3)-Za(1))-(Za(2)-Za(1))*(Ya(3)-Ya(1))
      ANy1=(Za(2)-Za(1))*(Xa(3)-Xa(1))-(Xa(2)-Xa(1))*(Za(3)-Za(1))
      ANz1=(Xa(2)-Xa(1))*(Za(3)-Za(1))-(Ya(2)-Ya(1))*(Xa(3)-Xa(1))
      ann(1)=(ANx1*ANx1+ANy1*ANy1+ANz1*ANz1)**0.5D0
 
      ANx2=(Ya(2)-Ya(1))*(Za(4)-Za(1))-(Za(2)-Za(1))*(Ya(4)-Ya(1))
      ANy2=(Za(2)-Za(1))*(Xa(4)-Xa(1))-(Xa(2)-Xa(1))*(Za(4)-Za(1))
      ANz2=(Xa(2)-Xa(1))*(Za(4)-Za(1))-(Ya(2)-Ya(1))*(Xa(4)-Xa(1))
      ann(2)=(ANx2*ANx2+ANy2*ANy2+ANz2*ANz2)**0.5D0
 
      IF(ann(1).GT.ann(2))THEN
      AN(1,NMS(1),K)=ANx1/ann(1)
      AN(2,NMS(1),K)=ANy1/ann(1)
      AN(3,NMS(1),K)=ANz1/ann(1)
      ELSE
      AN(1,NMS(1),K)=ANx2/ann(2)
      AN(2,NMS(1),K)=ANy2/ann(2)
      AN(3,NMS(1),K)=ANz2/ann(2)
      ann(1)=ann(2)
      ENDIF
 
      DO 531 I=2,NS
      AN(1,NMS(I),K)=AN(1,NMS(1),K)
      AN(2,NMS(I),K)=AN(2,NMS(1),K)
      AN(3,NMS(I),K)=AN(3,NMS(1),K)
 531  CONTINUE
 
c     write(6,*)'ann=',ann(1)
      DO 532 I=1,NY
      Xq=XYZ(1,NMY(I))
      Yq=XYZ(2,NMY(I))
      Zq=XYZ(3,NMY(I))
      YP(NMY(I),K)=ABS(
     &     (Xq-Xa(1))*AN(1,NMS(1),K)+(Yq-Ya(1))*AN(2,NMS(1),K)
     &                              +(Zq-Za(1))*AN(3,NMS(1),K))
c     write(6,*)'(Xq-Xa(1))=',(Xq-Xa(1)),' (Yq-Ya(1))=',(Yq-Ya(1)),
c    &' (Zq-Za(1))=',(Zq-Za(1))
 532  CONTINUE
 
c     write(6,*)'AN ',AN(1,NMS(1),K),AN(2,NMS(1),K),AN(3,NMS(1),K)
 
c     write(6,*)'NMS=',(NMS(I),I=1,NS)
c     write(6,*)'NMY=',(NMY(I),I=1,NY)
c     write(6,*)'YP=',(YP(NMS(I),K),I=1,NS)
c     write(6,*)'YP=',(YP(NMY(I),K),I=1,NY)
c     la normale n'est pas la meme pour tous les pts (SEG3)
      DO 533 I=1,NY
       AN(1,NMY(I),K)=AN(1,NMS(1),K)
       AN(2,NMY(I),K)=AN(2,NMS(1),K)
       AN(3,NMY(I),K)=AN(3,NMS(1),K)
 533  CONTINUE
c._._ FIN Calcul des normales aux CU27,PR21,TE15 et PY19
c._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._.
      ELSE
 
      write(6,*)'Type d elements non prevus dans FPU'
      write(6,*)'K=',K,' IDIM=',IDIM,' NS=',NS
 
      RETURN
      ENDIF
 
c     write(6,*)' FIN CALCUL NORMALES'
c._._ On Calcule la vitesse normale Un et la vitesse tangentielle Ut (NMY)
      DO 244 I=1,NP
      U0=0.D0
      U1=0.D0
      call initd(un,np,0.D0)
      DO 243 N=1,IDIM
      U0=U0+VS(I,N)*AN(N,I,K)
      U1=U1+(VS(I,N) - UN(I)*AN(N,I,K))**2.D0
 243  CONTINUE
      UN(I)=U0
      UT(I)=U1**0.5D0 + 1.D-10
 244  CONTINUE
c._._ On Calcule la vitesse normale Un et la vitesse tangentielle Ut (NMY)
 
c     write(6,*)'AN=',AN(1,NMS(1),K),AN(2,NMS(1),K)
c     write(6,*)'YP=',(YP(NMY(I),K),I=1,NY)
c     write(6,*)'UN=',(UN(NMY(I)),I=1,NY)
c     write(6,*)'UT=',(UT(NMY(I)),I=1,NY)
c     write(6,*)'VS= '
c     write(6,1002)(VS(I,1),I=1,NP)
c     write(6,1002)(VS(I,2),I=1,NP)
 
cc    CALL CALJBR(FN,GR,PG,XYZ,HR,PGSQ,RPG,
cc   &        NES,IDIM,NP,NPG,IAXI,AIRE,AJ,SGN)
c._._ On Calcule la matrice masse diagonale maillage appuye
cc    CALL INITD(DG,NPTA,0.D0)
cc    DO 221 I=1,NP
cc    AI=0.D0
cc    DO 220 LG=1,NPG
cc    PDR=PGSQ(LG)*DEUPI*RPG(LG)
cc    AI=AI+FN(I,LG)*PDR
c220  CONTINUE
cc    JD=IPADA(NUM(I,K))
cc    DG(JD)=DG(JD)+AI
c221  CONTINUE
c._._ On Calcule la matrice masse diagonale maillage appuye
 
      UETM=0.D0
      YPM =0.D0
      DO 225 II=1,NY
      I=NMY(II)
      IS=NMS(II)
      YP(I,K)=YP(I,K)+1.D-20
c     write(6,*)'YP ',YP(I,K)
 
      JI=IPADT(NUM(I,K))
c     JS=IPADS(NUM(IS,K))
      JR=(1-IKR)*(JI-1)+1
      JM=(1-IKM)*(JI-1)+1
c     YPLUS=RO(JR)*YP(I,K)*UETK(I,K)/AMU(JM)
      YPLUS=RO(JR)*YP(I,K)*UT(I)/AMU(JM)
      YPLUS=YPLUS+1.D-10
 
c     UET0=UETK(I,K)
      UET0=1.D-20
c     write(6,*)'M=0 YPLUS=',YPLUS,UT(I),AMU(JM),UET0
      w=0.5D0
 
      DO 30 M=1,NBIT
      IF(YPLUS.LE.5.D0)THEN
      UPLUS=YPLUS
      UETN=UT(I)/UPLUS
      YPLUS=w*YPLUS + (1.D0-w)*RO(JR)*YP(I,K)*UETN/AMU(JM)
      ELSE
c Reichardt
      UPLUS=1.D0/AKER*Log(1.D0+AKER*yplus)+
     & 7.8D0*(1.D0-exp(-yplus/11.D0)-(yplus/11.D0)*exp(-yplus/3.D0))
      UETN=UT(I)/UPLUS
      YPLUS=RO(JR)*YP(I,K)*UETN/AMU(JM)
      ENDIF
      resid=ABS(UETN-UET0)/UETN
      UET0=UETN
      IF(resid.LE.EPS)GO TO 31
 30   CONTINUE
      write(6,*)'Operateur FPU, NON CONVERGENCE en 40 iterations: eps=',
     & eps,' residu ->',resid
 31   CONTINUE
 
c     write(6,*)' Convergence en ',m,' iterations'
c     UETK(I,K)=UETN
c     VET(JS)=UETN
      UETM = (UETM*FLOAT(II-1)+UETN)/FLOAT(II)
      YPM  = (YPM *FLOAT(II-1)+YP(I,K))/FLOAT(II)
      AMUEF(I,K)=AMU(JM)
      IF(YPLUS.GT.5.D0)THEN
      AMUEF(I,K)=(RO(JR)*UETN*UETN*YP(I,K)/UT(I))+AMU(JM)
      ANUT(IPADT(NUM(I,K)))=AMUEF(I,K)
      ENDIF
 
 225  CONTINUE
c     write(6,*)' APRES 225 NY=',NY
 
      DO 226 I=1,NS
      JI=IPADT(NUM(I,K))
      JM=(1-IKM)*(JI-1)+1
      AMUEF(I,K)=AMU(JM)
      ANUT(IPADT(NUM(NMS(I),K)))=AMUEF(I,K)
      UET(IPADS(NUM(NMS(I),K)))=UETM
      YPS(IPADS(NUM(NMS(I),K)))=YPM
 226  CONTINUE
 
      IF(NBINC.EQ.3)THEN
      DO 227 I=1,NS
      JI=IPADT(NUM(I,K))
      JM=(1-IKM)*(JI-1)+1
      TKE(IPADS(NUM(NMS(I),K)),1)=UETM*UETM/SQCNU
      TKE(IPADS(NUM(NMS(I),K)),2)=UETM*UETM*UETM/(AKER*YPM)
 227  CONTINUE
      ENDIF
 
c       write(6,*)'NS=',NS
c       write(6,1001)(LKKE(I),I=1,NS)
c       write(6,1001)(NUM(I,K),I=1,NP)
c       write(6,*)'UET K=',K,(UETK(I,K),I=1,NP)
c       write(6,*)'AMUEF K=',K,(AMUEF(I,K),I=1,NP)
 
c._._ On deduit aussi Tau aux sommets
c._._ On deduit Nut aux points de Gauss
 
C.................................................................
 
C...... FIN Boucle sur les elements Indice K
c     write(6,*)'. FIN Boucle sur les elements Indice K ',K
  2      CONTINUE
 
c     write(6,*)'UET='
c     write(6,1001)(UET(i),i=1,npts)
 
 
c     write(6,*)' Nb d element traite KT=',KT
c     write(6,*)'RETOUR FPTAU'
      RETURN
 1002 FORMAT(10(1X,1PE11.4))
 1001 FORMAT(20(1X,I5))
      END
 
 
 
 

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