Sources Esope | Cast3M

kfnref
C KFNREF    SOURCE    GOUNAND   21/06/02    21:17:07     11022          
      SUBROUTINE KFNREF(MYLRF,MYPG,
     $     FNPG,DFNPG,
     $     IMPR,IRET)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : KFNREF
C PROJET      : Noyau linéaire NLIN
C DESCRIPTION : Calcul des fonctions de forme et de leurs
C               dérivées aux points de Gauss sur l'élément
C               de référence.
C               Le but avoué de ce sous-programme est d'effectuer
C               toutes les opérations de pré-traitement sur un élément
C               de référence donné.
C
C LANGAGE     : ESOPE
C AUTEUR      : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF)
C               mél : gounand@semt2.smts.cea.fr
C***********************************************************************
C APPELES          : CALPN, IVMAT, NI
C APPELE PAR       : KALPBG
C***********************************************************************
C ENTREES            : MYLRF, MYPG
C ENTREES/SORTIES    : -
C SORTIES            : FNPG, DFNPG
C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
C***********************************************************************
C VERSION    : v1,  26/07/99, version initiale
C HISTORIQUE : v1,  26/07/99, création
C HISTORIQUE : v1.1 16/09/99, remaniement pour que la base polynômiale
C                             puisse contenir des polynômes au lieu de
C                             monômes exclusivement...
C HISTORIQUE : v2, 10/05/00, modif. du segment ELREF
C HISTORIQUE : v2.1 02/10/03 ajout d'astuces
C HISTORIQUE :
C***********************************************************************
C Prière de PRENDRE LE TEMPS de compléter les commentaires
C en cas de modification de ce sous-programme afin de faciliter
C la maintenance !
C***********************************************************************
-INC PPARAM
-INC CCOPTIO
-INC CCREEL
-INC TNLIN      
*-INC SELREF
      POINTEUR MYLRF.ELREF
*-INC SPOGAU
      POINTEUR MYPG.POGAU
*-INC SMCHAEL
      POINTEUR FNPG.MCHEVA
      POINTEUR DFNPG.MCHEVA
-INC SMLENTI
      INTEGER JG
      POINTEUR IVTMP.MLENTI
-INC TMXMAT      
      POINTEUR PN.MXMAT
      POINTEUR PNM1.MXMAT
*
      INTEGER IMPR,IRET
*
      LOGICAL LLAHE,LCAST
      LOGICAL LCSTE,LLINE,LEGAL
      INTEGER NDFN
      REAL*8 DETPN
*
* Executable statements
*
      IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans kfnref'
      SEGACT MYLRF
*
      LCAST=(MYLRF.NOMLRF.EQ.'H1D1PY5'.OR.MYLRF.NOMLRF.EQ.'H1D2PY13'
     $     .OR.MYLRF.NOMLRF.EQ.'H1D2PR15'.OR.MYLRF.NOMLRF.EQ.'H1D2CU20')
*
      NDIML=MYLRF.ORDDER(/1)
      NBDDL=MYLRF.ORDDER(/2)
      LLAHE=(MYLRF.TYPEL.EQ.'LAGRANGE'.OR.MYLRF.TYPEL.EQ.'HERMITE')
* On fait un cas particulier aux éléments de Lagrange
* à un ddl : dans ce cas-la, la fonction de forme
* est constante et sa dérivée est nulle
      LCSTE=(MYLRF.TYPEL.EQ.'LAGRANGE'.AND.NBDDL.EQ.1)
* On fait un cas particulier aux éléments de Lagrange
* linéaire (simplex) : dans ce cas-la, les dérivées des fonctions de
* forme sont constantes
      LLINE=(MYLRF.TYPEL.EQ.'LAGRANGE'.AND.NBDDL.EQ.(NDIML+1))
      SEGDES MYLRF
*
* On repique les éléments dans SHAPE
*
      IF (LCAST) THEN
         CALL SH2FNR(MYLRF,MYPG,
     $        FNPG,DFNPG,
     $        IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
*
* Astuce foireuse
*
      ELSEIF (LCSTE) THEN
         NBLIG=1
         NBCOL=1
         N2LIG=1
         N2COL=1
         NBPOI=1
         NBELM=1
         SEGINI FNPG
         FNPG.WELCHE(1,1,1,1,1,1)=1.D0
         SEGDES FNPG
         NBLIG=1
         NBCOL=1
         N2LIG=1
         N2COL=NDIML
         NBPOI=1
         NBELM=1
         SEGINI DFNPG
         DO IDIML=1,NDIML
            DFNPG.WELCHE(1,1,1,IDIML,1,1)=0.D0
         ENDDO
         SEGDES DFNPG
* Cas des éléments de Lagrange et Hermite
      ELSEIF (LLAHE) THEN
*
*   Construisons la Matrice [PN] à l'aide de la base polynômiale :
*   et des coordonnées de noeuds d'approximation
**
         CALL CALPN(MYLRF,
     $        PN,
     $        IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
*
*   Inversons la Matrice [PN]
*
         SEGACT PN
         SEGINI,PNM1=PN
         NDFN=PN.XMAT(/1)
         JG=NDFN
         SEGINI IVTMP
         IIMPR=1
         CALL IVMAT(NDFN,PN.XMAT,
     $        IVTMP.LECT,
     $        PNM1.XMAT,DETPN,
     $        IIMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
         SEGSUP IVTMP
         SEGSUP PN
*
*   On peut maintenant calculer les valeurs des fonctions de forme
*   et leurs dérivées premières (par rapport aux coordonnées de l'espace
*   de référence) aux points de Gauss sur l'élément de référence
*
         CALL NI(MYLRF,MYPG,PNM1,
     $        FNPG,DFNPG,
     $        IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
         SEGSUP PNM1
      ELSE
         WRITE(IOIMP,*) 'Le type d''élément ',MYLRF.TYPEL
         WRITE(IOIMP,*) 'n''est pas reconnu.'
         GOTO 9999
      ENDIF
*
* Astuce foireuse (mais on fait quand même des
* vérifications)
*
      IF (LLINE) THEN
         SEGACT DFNPG*MOD
         NBPOGO=DFNPG.WELCHE(/5)
         DO IDIML=1,NDIML
            DO IDDL=1,NBDDL
               VALDF=DFNPG.WELCHE(1,IDDL,1,IDIML,1,1)
               DO IPOGO=2,NBPOGO
                  VALDF2=DFNPG.WELCHE(1,IDDL,1,IDIML,IPOGO,1)
                  LEGAL=(ABS(VALDF2-VALDF).LE.100.D0*XZPREC)
                  IF (.NOT.LEGAL) THEN
                     WRITE(IOIMP,*) 'Houston, on a un probleme'
*                     SEGPRT,DFNPG
                     GOTO 9999
                  ENDIF
               ENDDO
            ENDDO
         ENDDO
         NBLIG=1
         NBCOL=NBDDL
         N2LIG=1
         N2COL=NDIML
         NBPOI=1
         NBELM=1
         SEGADJ,DFNPG
         SEGDES DFNPG
      ENDIF
*
* Normal termination
*
      IRET=0
      RETURN
*
* Error handling
*
 9999 CONTINUE
      IRET=1
      WRITE(IOIMP,*) 'An error was detected in subroutine kfnref'
      RETURN
*
* End of subroutine KFNREF
*
      END