C NI SOURCE FANDEUR 22/05/02 21:15:26 11359 SUBROUTINE NI(MYLRF,MYPG,PNM1, $ FNPG,DFNPG, $ IMPR,IRET) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C*********************************************************************** C NOM : NI C PROJET : Noyau linéaire NLIN C DESCRIPTION : Calcul des valeurs des fonctions de forme et des C valeurs de leurs dérivées premières aux points de Gauss. C C On utilise la relation (produit) suivante : C C < Ni (point) > = < Pi (point) > [Pn]^{-1} C C avec (cf. CALPN) : C [PN] = ( P1(ksi1) ..... Pn(ksi1)) C ( ... ..... ... ) C ( P1(ksin) ..... Pn(ksin)) C n = nb. ddl sur l'élément (NDFN) C ksii = coords. du ieme noeud d'approximation C dans l'espace de référence (de dimension C NDIML) C Pi = ieme polynome d'interpolation sur C l'élément de référence. C Ni = ieme fonction nodale d'interpolation sur C l'élément de référence. C point= point quelconque sur l'élément de référence C (donc en particulier les points de Gauss) C LANGAGE : Esope C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF) C mél : gounand@semt2.smts.cea.fr C*********************************************************************** C APPELES : VALPOL C APPELE PAR : KFNREF C*********************************************************************** C ENTREES : MYLRF, MYPG, PNM1 C ENTREES/SORTIES : - C SORTIES : FNPG, DFNPG C CODE RETOUR (IRET) : = 0 si tout s'est bien passé C*********************************************************************** C VERSION : v1, 16/09/99, version initiale C HISTORIQUE : v1, 16/09/99, 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 SPOGAU POINTEUR MYPG.POGAU *-INC SELREF POINTEUR MYLRF.ELREF *-INC SPOLYNO POINTEUR MYBPOL.POLYNS POINTEUR MYPOLY.POLYNO -INC TMXMAT POINTEUR PNM1.MXMAT -INC SMLENTI POINTEUR ORDER1.MLENTI -INC SMLREEL POINTEUR VECTPI.MLREEL *-INC SMCHAEL INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM POINTEUR FNPG.MCHEVA POINTEUR DFNPG.MCHEVA * INTEGER IMPR,IRET * * Fonction Blas (produit scalaire) * REAL*8 DDOT EXTERNAL DDOT * INTEGER NBMONO INTEGER NDIML,NDPG,NBFN INTEGER INDIML,JNDIML,INDPG,INBFN * * Executable statements * IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans ni' * * Initialisations * SEGACT MYPG NDIML=MYPG.XCOPG(/1) NDPG=MYPG.XCOPG(/2) SEGACT MYLRF MYBPOL=MYLRF.MBPOLY SEGDES MYLRF SEGACT MYBPOL SEGACT MYBPOL.LIPOLY(*) NBFN=MYBPOL.LIPOLY(/1) SEGACT PNM1 JG=NDIML SEGINI ORDER1 JG=NBFN SEGINI VECTPI * * On calcule les valeurs des fonctions de forme aux points de Gauss * NBLIG=1 NBCOL=NBFN N2LIG=1 N2COL=1 NBPOI=NDPG NBELM=1 SEGINI FNPG DO 1 INDIML=1,NDIML ORDER1.LECT(INDIML)=0 1 CONTINUE DO 3 INDPG=1,NDPG * Calcul de < P (pg) > = < P1(pg) ... Pnbfn(pg) > où pg est le * INDPGieme point de Gauss DO 32 INBFN=1,NBFN MYPOLY=MYBPOL.LIPOLY(INBFN) NBMONO=MYPOLY.EXPMON(/2) CALL VALPOL(NDIML,NBMONO, $ MYPG.XCOPG(1,INDPG), $ MYPOLY.COEMON,MYPOLY.EXPMON, $ ORDER1.LECT, $ VECTPI.PROG(INBFN), $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 32 CONTINUE * On calcule : < N (pg) > = < P (pg) > [Pn]^{-1} DO 34 INBFN=1,NBFN FNPG.WELCHE(1,INBFN,1,1,INDPG,1)= $ DDOT(NBFN,VECTPI.PROG,1,PNM1.XMAT(1,INBFN),1) 34 CONTINUE 3 CONTINUE IF (IMPR.GT.3) THEN WRITE(IOIMP,*) 'Ordre de dérivation / coordonnée de réf. :' WRITE(IOIMP,4003) (ORDER1.LECT(INDIML),INDIML=1,NDIML) DO 5 INDPG=1,NDPG WRITE(IOIMP,*) 'Noeud de coordonnées :' WRITE(IOIMP,4004) (MYPG.XCOPG(INDIML,INDPG), $ INDIML=1,NDIML) WRITE(IOIMP,*) 'FNPG.WELCHE(nb.fns.forme) :' WRITE(IOIMP,4004) (FNPG.WELCHE(1,INBFN,1,1,INDPG,1), $ INBFN=1,NBFN) 5 CONTINUE ENDIF SEGDES FNPG * * On calcule les valeurs des dérivées premières des fonctions * de forme aux points de Gauss * NBLIG=1 NBCOL=NBFN N2LIG=1 N2COL=NDIML NBPOI=NDPG NBELM=1 SEGINI DFNPG IF (IRET.NE.0) GOTO 9999 DO 7 INDIML=1,NDIML DO 72 JNDIML=1,NDIML IF (JNDIML.EQ.INDIML) THEN ORDER1.LECT(JNDIML)=1 ELSE ORDER1.LECT(JNDIML)=0 ENDIF 72 CONTINUE DO 74 INDPG=1,NDPG DO 742 INBFN=1,NBFN * Calcul de < dP/dksi_indiml (pg) > où pg est le * INDPGieme point de Gauss MYPOLY=MYBPOL.LIPOLY(INBFN) NBMONO=MYPOLY.EXPMON(/2) CALL VALPOL(NDIML,NBMONO, $ MYPG.XCOPG(1,INDPG), $ MYPOLY.COEMON,MYPOLY.EXPMON, $ ORDER1.LECT, $ VECTPI.PROG(INBFN), $ IMPR,IRET) IF (IRET.NE.0) GOTO 9999 742 CONTINUE * On calcule : < dN/dksi_indiml (pg) > = < dP/dksi_indiml (pg) > [Pn]^{-1} DO 744 INBFN=1,NBFN DFNPG.WELCHE(1,INBFN,1,INDIML,INDPG,1)= $ DDOT(NBFN,VECTPI.PROG,1,PNM1.XMAT(1,INBFN),1) 744 CONTINUE 74 CONTINUE IF (IMPR.GT.3) THEN WRITE(IOIMP,*) $ 'Ordre de dérivation / coordonnée de réf. :' WRITE(IOIMP,4003) (ORDER1.LECT(JNDIML),JNDIML=1,NDIML) DO 76 INDPG=1,NDPG WRITE(IOIMP,*) 'Noeud de coordonnées :' WRITE(IOIMP,4004) (MYPG.XCOPG(JNDIML,INDPG), $ JNDIML=1,NDIML) WRITE(IOIMP,*) 'DFNPG.WELCHE(nb.fns.forme) :' WRITE(IOIMP,4004) $ (DFNPG.WELCHE(1,INBFN,1,INDIML,INDPG,1), $ INBFN=1,NBFN) 76 CONTINUE ENDIF 7 CONTINUE SEGDES DFNPG SEGSUP VECTPI SEGSUP ORDER1 SEGDES PNM1 SEGDES MYBPOL.LIPOLY(*) SEGDES MYBPOL SEGDES MYPG * * Normal termination * IRET=0 RETURN * * Format handling * 4003 FORMAT (2X,6(1X,I3)) 4004 FORMAT (2X,6(1X,1PE13.5)) * * Error handling * 9999 CONTINUE IRET=1 WRITE(IOIMP,*) 'An error was detected in subroutine ni' RETURN * * End of subroutine NI * END