Télécharger shapx.eso

Retour à la liste

shapx
C SHAPX     SOURCE    BP208322  18/04/12    21:15:12     9802           
c
      SUBROUTINE SHAPX(MELE,LV7,NBRMAX,NBENRJ,TLREEL,SHP,IRET)
C
C=======================================================================
C
C Entree: SHP = FONCTIONS DE FORME STANDARD ET LEUR DERIVEES
C         LV7 = FONCTIONS DE NIVEAU ET LEUR DERIVEES
C         LV7(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD)
C Sortie: SHP = TOUTES LES FONCTIONS DE FORME DE L'EF (std + enrichissement)
c               ET LEUR DERIVEES DANS LE REPERE GLOBAL
C
C
C=====PARTIE DECLARATIVE================================================
C
      implicit real*8(a-h,o-z)
      implicit integer(i-n)
 
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC SMLREEL
C
      DIMENSION   SHP(6,*)
c      PARAMETER  (NBRMAX=5)
      REAL*8      LV7(NBRMAX,2,6,*)
      INTEGER     TLREEL(NBRMAX,*)
c      DIMENSION   LV7(NBRMAX,2,6,*)
c      DIMENSION   TLREEL(NBRMAX,*)
      REAL*8      PSI,PHI,H,R,THETA
C
C
c      WRITE(*,*) '####    ENTREE DANS SHAPX   ####'
C
c
C=====AIGUILLAGE selon element==========================================
c
      if(MELE.eq.263) GOTO 263
      if(MELE.eq.264) GOTO 264
c     Erreur pour les elements non programmés
      call erreur (21)
      RETURN
C
CTY   if(MELE.eq.263) GOTO 263
CTY   Erreur pour les elements non programmés
CTY   call erreur (21)
CTY   RETURN
c
c
C=====XQ4R==============================================================
c
 263  continue
      NBNN = 4
c      KMAX =LV7(/1)
c
c
C++++ ENRICHISSEMENT H ++++++++++++++++++++++++++++++
      KENR = 1
 
C>>>> boucle sur les noeuds
      DO 200 I=1,NBNN
 
        MLREEL=TLREEL(KENR,I)
        if(MLREEL.eq.0) goto 200
 
C       on calcule les fonctions utiles
        PHI = LV7(KENR,2,1,I)
        H = SIGN(1.D0,PHI)
C
c       on modifie les fonctions de forme et leurs derivees
        II = (KENR*NBNN) + I
        DO IND=1,3
           SHP(IND,II) = H*SHP(IND,I)
        ENDDO
C
  200 CONTINUE
C<<<<<<<<<<<<<
C
C
C++++ ENRICHISSEMENT F1 et F2 ++++++++++++++++++++++++
C
c      DO 300 KENR=2,NBRMAX
      DO 300 KENR=2,NBENRJ
      IENR = 4*(KENR-2) + 2
 
C>>>>>>> boucle sur les noeuds
      DO 305 I=1,NBNN
C
        MLREEL=TLREEL(KENR,I)
c      write(*,*) 'KENR,I=',KENR,I,' MLREEL=',MLREEL
        if(MLREEL.eq.0) goto 305
 
C       on recupere les level set et leur derivees
        PSI  = LV7(KENR,1,1,I)
        PSIX = LV7(KENR,1,2,I)
        PSIY = LV7(KENR,1,3,I)
        PHI  = LV7(KENR,2,1,I)
        PHIX = LV7(KENR,2,2,I)
        PHIY = LV7(KENR,2,3,I)
 
C       on calcule les fonctions utiles
        H = SIGN(1.D0,PHI)
        R   = ( (PSI**2) + (PHI**2) )**0.5
        THETA = H * ((XPI/2.) - (ATAN(PSI/ABS(PHI))))
        SIN1T = SIN(THETA)
        COS1T = COS(THETA)
        SIN05T = SIN(THETA/2.)
        COS05T = COS(THETA/2.)
C
C       et leurs derivées
C       RX = ( (PSI*PSIX) + (PHI*PHIX) ) / (R**0.5)
C       RY = ( (PSI*PSIY) + (PHI*PHIY) ) / (R**0.5)
C       THETAX= ( ((-1.*ABS(PHI))*PSIX) + ((H*PSI)*PHIX) ) / (R**2)
C       THETAY= ( ((-1.*ABS(PHI))*PSIY) + ((H*PSI)*PHIY) ) / (R**2)
 
c       on calcule les fonctions de forme d ENRICHISSEMENT + derivees
 
C+++ Fb = (r**0.5)*sin(theta/2)
        II = (IENR*NBNN) + I
 
C       Fonction de forme Ni*Fb
        FB = (R**0.5) * SIN05T
        SHP(1,II) = SHP(1,I) * FB
C       Derivees (Ni*Fb),X  =  (Ni,X*Fb) + (Ni*Fb,X)
        SHP(2,II) = ( SHP(2,I) * FB ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( SIN05T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
     $     + ( COS05T * ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
C       Derivees (Ni*Fb),y  =  (Ni,y*Fb) + (Ni*Fb,y)
        SHP(3,II) = ( SHP(3,I) * FB ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( SIN05T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
     $     + ( COS05T * ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
 
 
C+++ Fc = (r**0.5)*cos(theta/2)
        II = ((IENR+1)*NBNN) + I
 
C       Fonction de forme Ni*Fc
        FC    = (R**0.5) * COS05T
        SHP(1,II) = SHP(1,I) * FC
C       Derivees (Ni*Fc),x  =  (Ni,x*Fc) + (Ni*Fc,x)
        SHP(2,II) = ( SHP(2,I) * FC ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( COS05T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
     $     - ( SIN05T * ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
C       Derivees (Ni*Fc),y  =  (Ni,y*Fc) + (Ni*Fc,y)
        SHP(3,II) = ( SHP(3,I) * FC) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( COS05T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
     $     - ( SIN05T * ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
 
 
C+++ Fd = (r**0.5)*sin(theta/2)*sin(theta)
        II = ((IENR+2)*NBNN) + I
 
C       Fonction de forme Ni*Fd
        FD    = (R**0.5) * SIN05T * SIN1T
        SHP(1,II)= SHP(1,I) * FD
C       Derivees (Ni*Fc),X  =  (Ni,X*Fc) + (Ni*Fc,X)
        SHP(2,II) = ( SHP(2,I) * FD )  +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( SIN05T*SIN1T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
     $     + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) *
     $                     ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
C       Derivees (Ni*Fc),Y  =  (Ni,Y*Fc) + (Ni*Fc,Y)
        SHP(3,II) = ( SHP(3,I) * FD ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( SIN05T*SIN1T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
     $     + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) *
     $                     ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
 
 
C+++ Fe = (r**0.5)*cos(theta/2)*sin(theta)
        II = ((IENR+3)*NBNN) + I
 
C       Fonction de forme Ni*Fe
        FE    = (R**0.5) * COS05T * SIN1T
        SHP(1,II)= SHP(1,I) * FE
C       Derivees (Ni*Fc),x  =  (Ni,x*Fe) + (Ni*Fe,x)
        SHP(2,II) = ( SHP(2,I) * FE ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( COS05T*SIN1T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
     $     - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) *
     $                     ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
C     Derivees (Ni*Fe),y  =  (Ni,y*Fe) + (Ni*Fe,y)
        SHP(3,II) = ( SHP(3,I) * FE ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( COS05T*SIN1T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
     $     - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) *
     $                     ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
 
C
C
  305 CONTINUE
C<<<<<<<<<<<<<
 
  300 CONTINUE
 
       GOTO 999
 
C_____________________________________________
C     FONCTIONS DE FORME TRIDIMENSIONNELLES
C
C=====XC8R==============================================================
c
 264  continue
C
      NBNN = 8
c      KMAX =LV7(/1)
c
c
C++++ ENRICHISSEMENT H ++++++++++++++++++++++++++++++
      KENR = 1
 
C>>>> boucle sur les noeuds
      DO 400 I=1,NBNN
 
        MLREEL=TLREEL(KENR,I)
        if(MLREEL.eq.0) goto 400
 
C       on calcule les fonctions utiles
        PHI = LV7(KENR,2,1,I)
        H = SIGN(1.D0,PHI)
C
c       on modifie les fonctions de forme et leurs derivees
        II = (KENR*NBNN) + I
        DO IND=1,4
           SHP(IND,II) = H*SHP(IND,I)
        ENDDO
C
  400 CONTINUE
C<<<<<<<<<<<<<
C
C
C++++ ENRICHISSEMENT F1 et F2 ++++++++++++++++++++++++
C
c      DO 500 KENR=2,NBRMAX
      DO 500 KENR=2,NBENRJ
      IENR = 4*(KENR-2) + 2
 
C>>>>>>> boucle sur les noeuds
      DO 505 I=1,NBNN
C
        MLREEL=TLREEL(KENR,I)
c      write(*,*) 'KENR,I=',KENR,I,' MLREEL=',MLREEL
        if(MLREEL.eq.0) goto 505
 
C       on recupere les level set et leur derivees
        PSI  = LV7(KENR,1,1,I)
        PSIX = LV7(KENR,1,2,I)
        PSIY = LV7(KENR,1,3,I)
        PSIZ = LV7(KENR,1,4,I)
        PHI  = LV7(KENR,2,1,I)
        PHIX = LV7(KENR,2,2,I)
        PHIY = LV7(KENR,2,3,I)
        PHIZ = LV7(KENR,2,4,I)
 
C       on calcule les fonctions utiles
        H = SIGN(1.D0,PHI)
        R   = ( (PSI**2) + (PHI**2) )**0.5
        THETA = H * ((XPI/2.) - (ATAN(PSI/ABS(PHI))))
        SIN1T = SIN(THETA)
        COS1T = COS(THETA)
        SIN05T = SIN(THETA/2.)
        COS05T = COS(THETA/2.)
C
C       et leurs derivées
C       RX = ( (PSI*PSIX) + (PHI*PHIX) ) / (R**0.5)
C       RY = ( (PSI*PSIY) + (PHI*PHIY) ) / (R**0.5)
C       THETAX= ( ((-1.*ABS(PHI))*PSIX) + ((H*PSI)*PHIX) ) / (R**2)
C       THETAY= ( ((-1.*ABS(PHI))*PSIY) + ((H*PSI)*PHIY) ) / (R**2)
 
c       on calcule les fonctions de forme d ENRICHISSEMENT + derivees
 
C+++ Fb = (r**0.5)*sin(theta/2)
        II = (IENR*NBNN) + I
 
C       Fonction de forme Ni*Fb
        FB = (R**0.5) * SIN05T
        SHP(1,II) = SHP(1,I) * FB
C       Derivees (Ni*Fb),X  =  (Ni,X*Fb) + (Ni*Fb,X)
        SHP(2,II) = ( SHP(2,I) * FB ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( SIN05T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
     $     + ( COS05T * ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
C       Derivees (Ni*Fb),y  =  (Ni,y*Fb) + (Ni*Fb,y)
        SHP(3,II) = ( SHP(3,I) * FB ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( SIN05T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
     $     + ( COS05T * ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
C       Derivees (Ni*Fb),z  =  (Ni,z*Fb) + (Ni*Fb,z)
        SHP(4,II) = ( SHP(4,I) * FB ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( SIN05T * ((COS1T*PSIZ)+(SIN1T*PHIZ)) )
     $     + ( COS05T * ((-1.*SIN1T*PSIZ)+(COS1T*PHIZ)) ) ) )
C
C+++ Fc = (r**0.5)*cos(theta/2)
        II = ((IENR+1)*NBNN) + I
 
C       Fonction de forme Ni*Fc
        FC    = (R**0.5) * COS05T
        SHP(1,II) = SHP(1,I) * FC
C       Derivees (Ni*Fc),x  =  (Ni,x*Fc) + (Ni*Fc,x)
        SHP(2,II) = ( SHP(2,I) * FC ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( COS05T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
     $     - ( SIN05T * ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
C       Derivees (Ni*Fc),y  =  (Ni,y*Fc) + (Ni*Fc,y)
        SHP(3,II) = ( SHP(3,I) * FC) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( COS05T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
     $     - ( SIN05T * ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
C       Derivees (Ni*Fc),z  =  (Ni,z*Fc) + (Ni*Fc,z)
        SHP(4,II) = ( SHP(4,I) * FC) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( COS05T * ((COS1T*PSIZ)+(SIN1T*PHIZ)) )
     $     - ( SIN05T * ((-1.*SIN1T*PSIZ)+(COS1T*PHIZ)) ) ) )
 
 
C+++ Fd = (r**0.5)*sin(theta/2)*sin(theta)
        II = ((IENR+2)*NBNN) + I
 
C       Fonction de forme Ni*Fd
        FD    = (R**0.5) * SIN05T * SIN1T
        SHP(1,II)= SHP(1,I) * FD
C       Derivees (Ni*Fd),X  =  (Ni,X*Fd) + (Ni*Fd,X)
        SHP(2,II) = ( SHP(2,I) * FD )  +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( SIN05T*SIN1T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
     $     + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) *
     $                     ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
C       Derivees (Ni*Fd),Y  =  (Ni,Y*Fd) + (Ni*Fd,Y)
        SHP(3,II) = ( SHP(3,I) * FD ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( SIN05T*SIN1T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
     $     + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) *
     $                     ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
C       Derivees (Ni*Fd),z  =  (Ni,z*Fd) + (Ni*Fd,z)
        SHP(4,II) = ( SHP(4,I) * FD ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( SIN05T*SIN1T * ((COS1T*PSIZ)+(SIN1T*PHIZ)) )
     $     + ( (COS05T*SIN1T + (2.*SIN05T*COS1T)) *
     $                     ((-1.*SIN1T*PSIZ)+(COS1T*PHIZ)) ) ) )
 
 
C+++ Fe = (r**0.5)*cos(theta/2)*sin(theta)
        II = ((IENR+3)*NBNN) + I
 
C       Fonction de forme Ni*Fe
        FE    = (R**0.5) * COS05T * SIN1T
        SHP(1,II)= SHP(1,I) * FE
C       Derivees (Ni*Fe),x  =  (Ni,x*Fe) + (Ni*Fe,x)
        SHP(2,II) = ( SHP(2,I) * FE ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( COS05T*SIN1T * ((COS1T*PSIX)+(SIN1T*PHIX)) )
     $     - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) *
     $                     ((-1.*SIN1T*PSIX)+(COS1T*PHIX)) ) ) )
C     Derivees (Ni*Fe),y  =  (Ni,y*Fe) + (Ni*Fe,y)
        SHP(3,II) = ( SHP(3,I) * FE ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( COS05T*SIN1T * ((COS1T*PSIY)+(SIN1T*PHIY)) )
     $     - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) *
     $                     ((-1.*SIN1T*PSIY)+(COS1T*PHIY)) ) ) )
C       Derivees (Ni*Fe),z  =  (Ni,z*Fe) + (Ni*Fe,z)
        SHP(4,II) = ( SHP(4,I) * FE ) +
     $  ( SHP(1,I)*0.5/(R**0.5) *
     $    (  ( COS05T*SIN1T * ((COS1T*PSIZ)+(SIN1T*PHIZ)) )
     $     - ( (SIN05T*SIN1T - (2.*COS05T*COS1T)) *
     $                     ((-1.*SIN1T*PSIZ)+(COS1T*PHIZ)) ) ) )
 
C
C
  505 CONTINUE
C<<<<<<<<<<<<<
 
  500 CONTINUE
 
 
C_____________________________________________
C     FIN DU PROGRAMME
  999 CONTINUE
C
 
      IRET = 1
c      WRITE(*,*) 'C EST FINI POUR SHAPX'
 
      RETURN
      END
 
 
 
 
 
 

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