C INELCU    SOURCE    GOUNAND   06/12/19    21:15:18     5612           
      SUBROUTINE INELCU(MYLRFS,IMPR,IRET)
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
C***********************************************************************
C NOM         : INELCU
C PROJET      : Noyau linéaire NLIN
C DESCRIPTION : Remplit le segment des éléments de référence
C               avec les éléments de référence de dimension 3,
C               de forme géométrique cubique.
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          : INILRF, FILRF, PROLRF, INILAG, PROBAP, GBAPCO
C APPELE PAR       : INLRFS
C***********************************************************************
C ENTREES            : -
C ENTREES/SORTIES    : MYLRFS
C SORTIES            : -
C CODE RETOUR (IRET) : = 0 si tout s'est bien passé
C***********************************************************************
C VERSION    : v1, 23/03/00, version initiale
C HISTORIQUE : v1, 23/03/00, création
C HISTORIQUE : v2, 10/05/00, modif. du segment ELREF
C HISTORIQUE :
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 CCOPTIO
CBEGININCLUDE SELREF
      SEGMENT ELREF
      CHARACTER*(LNNOM)  NOMLRF
      CHARACTER*(LNFORM) FORME 
      CHARACTER*(LNTYPL) TYPEL 
      CHARACTER*(LNESP)  ESPACE
      INTEGER DEGRE
      REAL*8 XCONOD(NDIMEL,NBNOD)
      INTEGER NPQUAF(NBDDL)
      INTEGER NUMCMP(NBDDL)
      INTEGER QUENOD(NBDDL)
      INTEGER ORDDER(NDIMEL,NBDDL)
      POINTEUR MBPOLY.POLYNS
      ENDSEGMENT
      SEGMENT ELREFS
      POINTEUR LISEL(0).ELREF
      ENDSEGMENT
CENDINCLUDE SELREF
      POINTEUR MYLRFS.ELREFS
      POINTEUR ELCOUR.ELREF
      POINTEUR ELPRO1.ELREF
      POINTEUR ELPRO2.ELREF
CBEGININCLUDE SPOLYNO
      SEGMENT POLYNO
      REAL*8  COEMON(NBMON)
      INTEGER EXPMON(NDIML,NBMON)
      ENDSEGMENT 
      SEGMENT POLYNS
      POINTEUR LIPOLY(NBPOLY).POLYNO
      ENDSEGMENT
CENDINCLUDE SPOLYNO
      POINTEUR MYBPOL.POLYNS
      POINTEUR MYBPO1.POLYNS
      POINTEUR MYBPO2.POLYNS
*
      INTEGER IMPR,IRET
* Elément de nom : L2D1CU4
      REAL*8 ZERO,UNS2,UN
      PARAMETER (ZERO=0.D0,UNS2=0.5D0,UN=1.D0)
*
      INTEGER INDDL
*     
* Executable statements
*
      IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans inelcu'
*
* Elément de nom : L2D0CU1
* Sur un cube : élément de Lagrange, fonction L2, approximation
*   nodale, espace de référence de dimension 3, 1 noeud, 1 degrés de
*   liberté, degré de l'approximation : 0
*       
* In INILRF : SEGINI ELCOUR
      CALL INILRF('L2D0CU1','CUBE','LAGRANGE','L2',
     $     3,1,1,0,
     $     ELCOUR,
     $     IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      CALL FILRF('L2D0QU1',MYLRFS,ELPRO1,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      CALL FILRF('L2D0SE1',MYLRFS,ELPRO2,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      CALL PROLRF(ELPRO1,ELPRO2,ELCOUR,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      ELCOUR.NPQUAF(1)=27
      ELCOUR.NUMCMP(1)=1
* Initialise la correspondance ddl-noeud+ord.der
      CALL INILAG(ELCOUR,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      SEGACT ELPRO1
      MYBPO1=ELPRO1.MBPOLY
      SEGDES ELPRO1
      SEGACT ELPRO2
      MYBPO2=ELPRO2.MBPOLY
      SEGDES ELPRO2
* Calcule la base polynômiale produit
      CALL PROBAP(MYBPO1,MYBPO2,MYBPOL,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      ELCOUR.MBPOLY=MYBPOL
      SEGDES ELCOUR
      MYLRFS.LISEL(**)=ELCOUR
*
* Elément de nom : L2D1CU4
* Sur un cube : élément de Lagrange, fonction L2, approximation
*   nodale, espace de référence de dimension 3, 4 noeuds, 4 degrés de
*   liberté, degré de l'approximation : 1
*       
* In INILRF : SEGINI ELCOUR
      CALL INILRF('L2D1CU4','CUBE','LAGRANGE','L2',
     $     3,4,4,1,
     $     ELCOUR,
     $     IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      ELCOUR.XCONOD(1,1)= UNS2
      ELCOUR.XCONOD(2,1)= UNS2
      ELCOUR.XCONOD(3,1)=-UNS2
      ELCOUR.XCONOD(1,2)=-UNS2
      ELCOUR.XCONOD(2,2)= ZERO 
      ELCOUR.XCONOD(3,2)=-UNS2
      ELCOUR.XCONOD(1,3)= ZERO
      ELCOUR.XCONOD(2,3)=-UNS2
      ELCOUR.XCONOD(3,3)=-UNS2
      ELCOUR.XCONOD(1,4)= ZERO
      ELCOUR.XCONOD(2,4)=-UNS2
      ELCOUR.XCONOD(3,4)= UNS2
      ELCOUR.NPQUAF(1)=27
      ELCOUR.NUMCMP(1)=1
      ELCOUR.NPQUAF(2)=27
      ELCOUR.NUMCMP(2)=2
      ELCOUR.NPQUAF(3)=27
      ELCOUR.NUMCMP(3)=3
      ELCOUR.NPQUAF(4)=27
      ELCOUR.NUMCMP(4)=4
* Initialise la correspondance ddl-noeud+ord.der
      CALL INILAG(ELCOUR,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
* Génère une base polynômiale complète (dimension 3, degré 1)
      CALL GBAPCO(3,1,MYBPOL,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      ELCOUR.MBPOLY=MYBPOL
      SEGDES ELCOUR
      MYLRFS.LISEL(**)=ELCOUR
*
* Elément de nom : CRD1CU6
* Sur un cube : élément de Lagrange, fonction continue au centre
* des faces, approximation
*   nodale, espace de référence de dimension 3, 6 noeuds, 6 degrés de
*   liberté, degré de l'approximation : 1
*       
* In INILRF : SEGINI ELCOUR
      CALL INILRF('CRD1CU6','CUBE','LAGRANGE','HFAC',
     $     3,6,6,1,
     $     ELCOUR,
     $     IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      ELCOUR.XCONOD(1,1)=ZERO
      ELCOUR.XCONOD(2,1)=-UN
      ELCOUR.XCONOD(3,1)=ZERO
      ELCOUR.XCONOD(1,2)=UN
      ELCOUR.XCONOD(2,2)=ZERO
      ELCOUR.XCONOD(3,2)=ZERO
      ELCOUR.XCONOD(1,3)=ZERO
      ELCOUR.XCONOD(2,3)=UN
      ELCOUR.XCONOD(3,3)=ZERO
      ELCOUR.XCONOD(1,4)=-UN
      ELCOUR.XCONOD(2,4)=ZERO
      ELCOUR.XCONOD(3,4)=ZERO
      ELCOUR.XCONOD(1,5)=ZERO
      ELCOUR.XCONOD(2,5)=ZERO
      ELCOUR.XCONOD(3,5)=-UN
      ELCOUR.XCONOD(1,6)=ZERO
      ELCOUR.XCONOD(2,6)=ZERO
      ELCOUR.XCONOD(3,6)=+UN
* Les d.d.l. sont aux noeuds 21 à 26...
      DO INDDL=1,6
         ELCOUR.NPQUAF(INDDL)=INDDL+20
         ELCOUR.NUMCMP(INDDL)=1
      ENDDO
* Initialise la correspondance ddl-noeud+ord.der
      CALL INILAG(ELCOUR,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
* Génère une base polynômiale complète (dimension 3, degré 1)
      CALL GBAPCO(3,1,MYBPOL,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
* On rajoute les polynômes spécifiques à crouzeix-raviart quadrilatère
      CALL GPOCRQ(3,MYBPOL,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      ELCOUR.MBPOLY=MYBPOL
      SEGDES ELCOUR
      MYLRFS.LISEL(**)=ELCOUR
*
* Elément de nom : H1D1CU8
* Sur un cube : élément de Lagrange, fonction H1, approximation
*   nodale, espace de référence de dimension 3, 8 noeuds, 8 degrés de
*   liberté, degré de l'approximation : 1
*       
* In INILRF : SEGINI ELCOUR
      CALL INILRF('H1D1CU8','CUBE','LAGRANGE','H1',
     $     3,8,8,1,
     $     ELCOUR,
     $     IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      CALL FILRF('H1D1QU4',MYLRFS,ELPRO1,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      CALL FILRF('H1D1SE2',MYLRFS,ELPRO2,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      CALL PROLRF(ELPRO1,ELPRO2,ELCOUR,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
* Les d.d.l. sont aux noeuds 1,3,5,7...
      DO 201 INDDL=1,4
         ELCOUR.NPQUAF(INDDL)=(2*INDDL)-1
         ELCOUR.NUMCMP(INDDL)=1
 201  CONTINUE
* ... et 13,15,17,19.
      DO 203 INDDL=5,8
         ELCOUR.NPQUAF(INDDL)=(2*(INDDL-5))+13
         ELCOUR.NUMCMP(INDDL)=1
 203  CONTINUE
* Initialise la correspondance ddl-noeud+ord.der
      CALL INILAG(ELCOUR,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      SEGACT ELPRO1
      MYBPO1=ELPRO1.MBPOLY
      SEGDES ELPRO1
      SEGACT ELPRO2
      MYBPO2=ELPRO2.MBPOLY
      SEGDES ELPRO2
* Calcule la base polynômiale produit
      CALL PROBAP(MYBPO1,MYBPO2,MYBPOL,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      ELCOUR.MBPOLY=MYBPOL
      SEGDES ELCOUR
      MYLRFS.LISEL(**)=ELCOUR
*
* Elément de nom : H1D2CU20
* Sur un cube : élément de Lagrange, fonction H1, approximation
*   nodale, espace de référence de dimension 3, 20 noeuds, 20 degrés de
*   liberté, degré de l'approximation : 2
*       
* In INILRF : SEGINI ELCOUR
      CALL INILRF('H1D2CU20','CUBE','LAGRANGE','H1',
     $     3,20,20,2,
     $     ELCOUR,
     $     IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
C Inutile
C      ELCOUR.XCONOD(1,1)=ZERO
C      ELCOUR.XCONOD(2,1)=ZERO
C      ELCOUR.XCONOD(3,1)=ZERO
* Les d.d.l sont aux noeuds 1,...,20
      DO IDDL=1,20
         ELCOUR.NPQUAF(IDDL)=IDDL
         ELCOUR.NUMCMP(IDDL)=1
      ENDDO
* Initialise la correspondance ddl-noeud+ord.der
      CALL INILAG(ELCOUR,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
* Pas de base polynômiale (on recopie l'élément de castem)
      ELCOUR.MBPOLY=0
      SEGDES ELCOUR
      MYLRFS.LISEL(**)=ELCOUR
*
* Elément de nom : H1D2CU27
* Sur un cube : élément de Lagrange, fonction H1, approximation
*   nodale, espace de référence de dimension 3, 27 noeuds, 27 degrés de
*   liberté, degré de l'approximation : 2
*       
* In INILRF : SEGINI ELCOUR
      CALL INILRF('H1D2CU27','CUBE','LAGRANGE','H1',
     $     3,27,27,2,
     $     ELCOUR,
     $     IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      CALL FILRF('H1D2QU9',MYLRFS,ELPRO1,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      CALL FILRF('H1D2SE3',MYLRFS,ELPRO2,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      CALL PROLRF(ELPRO1,ELPRO2,ELCOUR,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
* Les ieme d.d.l sont aux noeuds j
      DO IDDL=1,8
         ELCOUR.NPQUAF(IDDL)=IDDL
      ENDDO
      ELCOUR.NPQUAF( 9)=25
      ELCOUR.NPQUAF(10)=9
      ELCOUR.NPQUAF(11)=21
      ELCOUR.NPQUAF(12)=10
      ELCOUR.NPQUAF(13)=22
      ELCOUR.NPQUAF(14)=11
      ELCOUR.NPQUAF(15)=23
      ELCOUR.NPQUAF(16)=12
      ELCOUR.NPQUAF(17)=24
      ELCOUR.NPQUAF(18)=27
      DO IDDL=19,26
         ELCOUR.NPQUAF(IDDL)=IDDL-6
      ENDDO
      ELCOUR.NPQUAF(27)=26
      DO IDDL=1,27
         ELCOUR.NUMCMP(IDDL)=1
      ENDDO
* Initialise la correspondance ddl-noeud+ord.der
      CALL INILAG(ELCOUR,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      SEGACT ELPRO1
      MYBPO1=ELPRO1.MBPOLY
      SEGDES ELPRO1
      SEGACT ELPRO2
      MYBPO2=ELPRO2.MBPOLY
      SEGDES ELPRO2
* Calcule la base polynômiale produit
      CALL PROBAP(MYBPO1,MYBPO2,MYBPOL,IMPR,IRET)
      IF (IRET.NE.0) GOTO 9999
      ELCOUR.MBPOLY=MYBPOL
      SEGDES ELCOUR
      MYLRFS.LISEL(**)=ELCOUR
*
* Normal termination
*
      IRET=0
      RETURN
*
* Format handling
*
*
* Error handling
*
 9999 CONTINUE
      IRET=1
      WRITE(IOIMP,*) 'An error was detected in subroutine inelcu'
      RETURN
*
* End of subroutine INELCU
*
      END