C INELCU SOURCE GOUNAND 21/06/02 21:16:26 11022 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 PPARAM -INC CCOPTIO -INC TNLIN *-INC SELREF POINTEUR MYLRFS.ELREFS POINTEUR ELCOUR.ELREF POINTEUR ELPRO1.ELREF POINTEUR ELPRO2.ELREF *-INC 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