C GEOLI2 SOURCE GOUNAND 26/01/09 21:15:23 12441 SUBROUTINE GEOLI2(IESREF,NDPOGO,NDELEM,JMAJAC, $ JMIJAC,JDTJAC, $ IMPR,IRET) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER (I-N) C*********************************************************************** C NOM : GEOLI2 C PROJET : Noyau linéaire NLIN C DESCRIPTION : Calcul de l'inverse et du déterminant d'une matrice C jacobienne dans le cas où celle-ci est carrée. C Ceci est effectué pour chaque point de Gauss d'un C élément. C On teste également si le déterminant s'annule ou change C de signe sur l'élément auxquels cas on génère une C erreur. C C LANGAGE : Fortran 77 (sauf E/S) C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF) C mél : gounand@semt2.smts.cea.fr C*********************************************************************** C APPELES : - C APPELE PAR : GEOLIN C*********************************************************************** C ENTREES : * IESREF (type entier) : dimension de l'espace de C référence. C * NDPOGO (type entier) : nombre de points C d'intégration. C * NDELEM (type entier) : nombre d'éléments du C maillage élémentaire courant. C * JMAJAC (type MCHEVA) : valeurs de la matrice C jacobienne aux points de Gauss sur le maillage C élémentaire courant. C ENTREES/SORTIES : * JMIJAC (type MCHEVA) : valeurs de l'inverse de C la matrice jacobienne aux points de Gauss sur C le maillage élémentaire courant. C JMIJAC n'existe que si dim.esp.réf=dim.esp.réel C * JDTJAC (type MCHEVA) : valeurs du déterminant C de la matrice jacobienne aux points de Gauss C sur le maillage élémentaire courant. C Si la matrice jacobienne n'est pas carrée, on C a calculé sqrt(det(trans(J).J)) C SORTIES : - C CODE RETOUR (IRET) : = 0 si tout s'est bien passé C*********************************************************************** C VERSION : v1, 10/08/99, version initiale C HISTORIQUE : v1, 10/08/99, création 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 CCREEL -INC PPARAM -INC CCOPTIO INTEGER IESREF,NDPOGO,NDELEM REAL*8 JMAJAC(IESREF,IESREF,NDPOGO,NDELEM) REAL*8 JMIJAC(IESREF,IESREF,NDPOGO,NDELEM) REAL*8 JDTJAC(NDPOGO,NDELEM) * INTEGER IMPR,IRET * REAL*8 ZERO,UN PARAMETER (ZERO=0.D0,UN=1.D0) * REAL*8 DET,INVDET INTEGER IELEM,IPOGO,IREFER,JREFER LOGICAL LSIGN,LDETZ,LCHSIG * * Executable statements * * La ligne suivante est inutile mais le compilateur du jour bugge si * on ne la mat pas LSIGN=.FALSE. LDETZ=.FALSE. LCHSIG=.FALSE. IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans geoli2' IF (IESREF.EQ.1) THEN DO 1 IELEM=1,NDELEM DO 12 IPOGO=1,NDPOGO DET=JMAJAC(1,1,IPOGO,IELEM) JDTJAC(IPOGO,IELEM)=DET * IF (DET.NE.ZERO) THEN IF (ABS(DET).GT.(SQRT(XPETIT))) THEN IF (IPOGO.EQ.1) THEN LSIGN=(DET.GT.ZERO) ELSE IF (LSIGN.NEQV.(DET.GT.ZERO)) LCHSIG=.TRUE. ENDIF INVDET=UN/DET JMIJAC(1,1,IPOGO,IELEM)=INVDET ELSE LDETZ=.TRUE. ENDIF 12 CONTINUE 1 CONTINUE ELSEIF (IESREF.EQ.2) THEN DO 3 IELEM=1,NDELEM DO 32 IPOGO=1,NDPOGO DET=(JMAJAC(1,1,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM)) $ - (JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(1,2,IPOGO,IELEM)) JDTJAC(IPOGO,IELEM)=DET * IF (DET.NE.ZERO) THEN IF (ABS(DET).GT.(SQRT(XPETIT))) THEN IF (IPOGO.EQ.1) THEN LSIGN=(DET.GT.ZERO) ELSE IF (LSIGN.NEQV.(DET.GT.ZERO)) LCHSIG=.TRUE. ENDIF INVDET=UN/DET JMIJAC(1,1,IPOGO,IELEM)= $ (+JMAJAC(2,2,IPOGO,IELEM))*INVDET JMIJAC(2,1,IPOGO,IELEM)= $ (-JMAJAC(2,1,IPOGO,IELEM))*INVDET JMIJAC(1,2,IPOGO,IELEM)= $ (-JMAJAC(1,2,IPOGO,IELEM))*INVDET JMIJAC(2,2,IPOGO,IELEM)= $ (+JMAJAC(1,1,IPOGO,IELEM))*INVDET ELSE LDETZ=.TRUE. ENDIF 32 CONTINUE 3 CONTINUE ELSEIF (IESREF.EQ.3) THEN DO 5 IELEM=1,NDELEM DO 52 IPOGO=1,NDPOGO JMIJAC(1,1,IPOGO,IELEM)= $ (JMAJAC(2,2,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM)) $ -(JMAJAC(3,2,IPOGO,IELEM)*JMAJAC(2,3,IPOGO,IELEM)) JMIJAC(2,1,IPOGO,IELEM)= $ (JMAJAC(3,1,IPOGO,IELEM)*JMAJAC(2,3,IPOGO,IELEM)) $ -(JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM)) JMIJAC(3,1,IPOGO,IELEM)= $ (JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(3,2,IPOGO,IELEM)) $ -(JMAJAC(3,1,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM)) JMIJAC(1,2,IPOGO,IELEM)= $ (JMAJAC(1,3,IPOGO,IELEM)*JMAJAC(3,2,IPOGO,IELEM)) $ -(JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM)) JMIJAC(2,2,IPOGO,IELEM)= $ (JMAJAC(1,1,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM)) $ -(JMAJAC(1,3,IPOGO,IELEM)*JMAJAC(3,1,IPOGO,IELEM)) JMIJAC(3,2,IPOGO,IELEM)= $ (JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(3,1,IPOGO,IELEM)) $ -(JMAJAC(3,2,IPOGO,IELEM)*JMAJAC(1,1,IPOGO,IELEM)) JMIJAC(1,3,IPOGO,IELEM)= $ (JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(2,3,IPOGO,IELEM)) $ -(JMAJAC(1,3,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM)) JMIJAC(2,3,IPOGO,IELEM)= $ (JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(1,3,IPOGO,IELEM)) $ -(JMAJAC(2,3,IPOGO,IELEM)*JMAJAC(1,1,IPOGO,IELEM)) JMIJAC(3,3,IPOGO,IELEM)= $ (JMAJAC(1,1,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM)) $ -(JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(2,1,IPOGO,IELEM)) DET=(JMAJAC(1,1,IPOGO,IELEM)*JMIJAC(1,1,IPOGO,IELEM)) $ +(JMAJAC(1,2,IPOGO,IELEM)*JMIJAC(2,1,IPOGO,IELEM)) $ +(JMAJAC(1,3,IPOGO,IELEM)*JMIJAC(3,1,IPOGO,IELEM)) JDTJAC(IPOGO,IELEM)=DET * IF (DET.NE.ZERO) THEN IF (ABS(DET).GT.(SQRT(XPETIT))) THEN IF (IPOGO.EQ.1) THEN LSIGN=(DET.GT.ZERO) ELSE IF (LSIGN.NEQV.(DET.GT.ZERO)) LCHSIG=.TRUE. ENDIF INVDET=UN/DET DO 522 JREFER=1,IESREF DO 5222 IREFER=1,IESREF JMIJAC(IREFER,JREFER,IPOGO,IELEM)= $ JMIJAC(IREFER,JREFER,IPOGO,IELEM) $ * INVDET 5222 CONTINUE 522 CONTINUE ELSE LDETZ=.TRUE. ENDIF 52 CONTINUE 5 CONTINUE ELSE WRITE(IOIMP,*) 'Je ne sais pas calculer l''inverse d''une' WRITE(IOIMP,*) 'matrice jacobienne de dimension ' WRITE(IOIMP,*) 'IESREF=',IESREF ENDIF IF (LDETZ) GOTO 9998 IF (LCHSIG) GOTO 9997 * * Normal termination * IRET=0 RETURN * * Format handling * * * Error handling * 9997 CONTINUE * IF (LERJ) THEN IRET=-666 * CALL ECRENT(IRET) RETURN * ENDIF * WRITE(IOIMP,*) 'Le déterminant de la matrice jacobienne' * WRITE(IOIMP,*) 'change de signe sur l''élément.' * WRITE(IOIMP,*) 'IELEM=',IELEM,' IPOGO=',IPOGO * GOTO 9999 9998 CONTINUE * IF (LERJ) THEN IRET=-667 * CALL ECRENT(IRET) RETURN * ENDIF * WRITE(IOIMP,*) 'Déterminant de la matrice jacobienne nul' * WRITE(IOIMP,*) 'IELEM=',IELEM,' IPOGO=',IPOGO * GOTO 9999 9999 CONTINUE IRET=1 WRITE(IOIMP,*) 'An error was detected in subroutine geoli2' RETURN * * End of subroutine GEOLI2 * END