Sources Esope | Cast3M

kfnrff
C KFNRFF    SOURCE    GOUNAND   21/06/02    21:17:08     11022          
      SUBROUTINE KFNRFF(LRFVOL,JXCOPG,
     $     FNPG,DFNPG,
     $     IMPR,IRET)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : KFNRFF
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            : LRFVOL, 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 LRFVOL.ELREF
*-INC SMCHAEL
      POINTEUR JXCOPG.MCHEVA
      POINTEUR FNPG.MCHEVA
      POINTEUR DFNPG.MCHEVA
-INC SMLENTI
      INTEGER JG
      POINTEUR IVTMP.MLENTI
-INC TMXMAT
      POINTEUR PN.MXMAT
      POINTEUR PNM1.MXMAT
      POINTEUR ID.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 kfnrff'
      SEGACT LRFVOL
*
      LCAST=(LRFVOL.NOMLRF.EQ.'H1D1PY5'.OR.LRFVOL.NOMLRF.EQ.'H1D2PY13'
     $     .OR.LRFVOL.NOMLRF.EQ.'H1D2PR15'
     $     .OR.LRFVOL.NOMLRF.EQ.'H1D2CU20')
*
      NDIML=LRFVOL.ORDDER(/1)
      NBDDL=LRFVOL.ORDDER(/2)
      LLAHE=(LRFVOL.TYPEL.EQ.'LAGRANGE'.OR.LRFVOL.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=(LRFVOL.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=(LRFVOL.TYPEL.EQ.'LAGRANGE'.AND.NBDDL.EQ.(NDIML+1))
      SEGDES LRFVOL
*
* On repique les éléments dans SHAPE
*
      IF (LCAST) THEN
         CALL SH2FNF(LRFVOL,JXCOPG,
     $        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(LRFVOL,
     $        PN,
     $        IMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
*
*   Inversons la Matrice [PN]
*
C       SEGACT PN
C       JLIG=PN.MAT(/1)
C       JCOL=PN.MAT(/2)
C       xsom=0
C       do icol=1,jcol
C          do ilig=1,jlig
C             xsom=xsom+ABS(PN.MAT(ilig,icol))
C          enddo
C       enddo
C       write(ioimp,*) 'xsom=',xsom
 
 
         SEGACT PN
         SEGINI,PNM1=PN
         NDFN=PN.XMAT(/1)
         JG=NDFN
         SEGINI IVTMP
         IIMPR=1
*         IIMPR=5
         CALL IVMAT(NDFN,PN.XMAT,
     $        IVTMP.LECT,
     $        PNM1.XMAT,DETPN,
     $        IIMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
C* Vérif à la con
C         JLIG=NDFN
C         JCOL=NDFN
C         SEGINI ID
C         DO i = 1,NDFN
C            DO k = 1,NDFN
C               val = 0.D0
C               DO j = 1,NDFN
C                  val = val + (PN.XMAT(i,j)*PNM1.XMAT(j,k))
C               ENDDO
C               ID.MAT(i,k)=val
C            enddo
C         enddo
C         segprt,ID
C         xsom=0.D0
C         do icol=1,jcol
C            do ilig=1,jlig
C               xsom=xsom+ABS(ID.MAT(ilig,icol))
C            enddo
C         enddo
C         write(ioimp,*) 'xsom=',xsom
C         write(ioimp,*) 'diagonale'
C         WRITE(IOIMP,*) (ID.MAT(I,I),I=1,NDFN)
C         xsom=0
C         do icol=1,jcol
C            do ilig=1,jlig
C               xsom=xsom+ABS(PNM1.MAT(ilig,icol))
C            enddo
C         enddo
C         write(ioimp,*) 'xsom=',xsom
         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
*
         IIMPR=IMPR
         CALL NIF(LRFVOL,JXCOPG,PNM1,
     $        FNPG,DFNPG,
     $        IIMPR,IRET)
         IF (IRET.NE.0) GOTO 9999
         SEGSUP PNM1
      ELSE
         WRITE(IOIMP,*) 'Le type d''élément ',LRFVOL.TYPEL
         WRITE(IOIMP,*) 'n''est pas reconnu.'
         GOTO 9999
      ENDIF
*
* Astuce foireuse (mais on fait quand même des
* vérifications) On ne la fait pas pour les intégrations de surface
*
C      IF (LLINE) THEN
C         SEGACT DFNPG*MOD
C         NBPOGO=DFNPG.WELCHE(/5)
C         NBELFV=DFNPG.WELCHE(/6)
C         DO IDIML=1,NDIML
C            DO IDDL=1,NBDDL
C               VALDF=DFNPG.WELCHE(1,IDDL,1,IDIML,1,1)
C               DO IBELFV=2,NBELFV
C                  DO IPOGO=2,NBPOGO
C                     VALDF2=DFNPG.WELCHE(1,IDDL,1,IDIML,IPOGO,1)
C                     LEGAL=(ABS(VALDF2-VALDF).LE.100.D0*XZPREC)
C                     IF (.NOT.LEGAL) THEN
C                        WRITE(IOIMP,*) 'Houston, on a un probleme'
C*                     SEGPRT,DFNPG
C                        GOTO 9999
C                     ENDIF
C                  ENDDO
C               ENDDO
C            ENDDO
C         ENDDO
C         NBLIG=1
C         NBCOL=NBDDL
C         N2LIG=1
C         N2COL=NDIML
C         NBPOI=1
C         NBELM=1
C         SEGADJ,DFNPG
C         SEGDES DFNPG
C      ENDIF
*
* Normal termination
*
      IRET=0
      RETURN
*
* Error handling
*
 9999 CONTINUE
      IRET=1
      WRITE(IOIMP,*) 'An error was detected in subroutine kfnrff'
      RETURN
*
* End of subroutine KFNRFF
*
      END