1 : C GEOLI2 SOURCE GOUNAND 05/12/21 21:27:00 5281 2 : SUBROUTINE GEOLI2(IESREF,NDPOGO,NDELEM,JMAJAC, 3 : $ JMIJAC,JDTJAC,LERJ, 4 : $ IMPR,IRET) 5 : IMPLICIT REAL*8 (A-H,O-Z) 6 : IMPLICIT INTEGER (I-N) 7 : C*********************************************************************** 8 : C NOM : GEOLI2 9 : C PROJET : Noyau linéaire NLIN 10 : C DESCRIPTION : Calcul de l'inverse et du déterminant d'une matrice 11 : C jacobienne dans le cas où celle-ci est carrée. 12 : C Ceci est effectué pour chaque point de Gauss d'un 13 : C élément. 14 : C On teste également si le déterminant s'annule ou change 15 : C de signe sur l'élément auxquels cas on génère une 16 : C erreur. 17 : C 18 : C LANGAGE : Fortran 77 (sauf E/S) 19 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF) 20 : C mél : gounand@semt2.smts.cea.fr 21 : C*********************************************************************** 22 : C APPELES : - 23 : C APPELE PAR : GEOLIN 24 : C*********************************************************************** 25 : C ENTREES : * IESREF (type entier) : dimension de l'espace de 26 : C référence. 27 : C * NDPOGO (type entier) : nombre de points 28 : C d'intégration. 29 : C * NDELEM (type entier) : nombre d'éléments du 30 : C maillage élémentaire courant. 31 : C * JMAJAC (type MCHEVA) : valeurs de la matrice 32 : C jacobienne aux points de Gauss sur le maillage 33 : C élémentaire courant. 34 : C ENTREES/SORTIES : * JMIJAC (type MCHEVA) : valeurs de l'inverse de 35 : C la matrice jacobienne aux points de Gauss sur 36 : C le maillage élémentaire courant. 37 : C JMIJAC n'existe que si dim.esp.réf=dim.esp.réel 38 : C * JDTJAC (type MCHEVA) : valeurs du déterminant 39 : C de la matrice jacobienne aux points de Gauss 40 : C sur le maillage élémentaire courant. 41 : C Si la matrice jacobienne n'est pas carrée, on 42 : C a calculé sqrt(det(trans(J).J)) 43 : C SORTIES : - 44 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé 45 : C*********************************************************************** 46 : C VERSION : v1, 10/08/99, version initiale 47 : C HISTORIQUE : v1, 10/08/99, création 48 : C HISTORIQUE : 49 : C HISTORIQUE : 50 : C*********************************************************************** 51 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 52 : C en cas de modification de ce sous-programme afin de faciliter 53 : C la maintenance ! 54 : C*********************************************************************** 55 : -INC CCREEL 56 : -INC CCOPTIO 57 : INTEGER IESREF,NDPOGO,NDELEM 58 : REAL*8 JMAJAC(IESREF,IESREF,NDPOGO,NDELEM) 59 : REAL*8 JMIJAC(IESREF,IESREF,NDPOGO,NDELEM) 60 : REAL*8 JDTJAC(NDPOGO,NDELEM) 61 : * 62 : INTEGER IMPR,IRET 63 : * 64 : REAL*8 ZERO,UN 65 : PARAMETER (ZERO=0.D0,UN=1.D0) 66 : * 67 : REAL*8 DET,INVDET 68 : INTEGER IELEM,IPOGO,IREFER,JREFER 69 : LOGICAL LSIGN,LERJ 70 : * 71 : * Executable statements 72 : * 73 : IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans geoli2' 74 : IF (IESREF.EQ.1) THEN 75 : DO 1 IELEM=1,NDELEM 76 : DO 12 IPOGO=1,NDPOGO 77 : DET=JMAJAC(1,1,IPOGO,IELEM) 78 : IF (DET.NE.ZERO) THEN 79 : INVDET=UN/DET 80 : JDTJAC(IPOGO,IELEM)=DET 81 : JMIJAC(1,1,IPOGO,IELEM)=INVDET 82 : ELSE 83 : GOTO 9998 84 : ENDIF 85 : 12 CONTINUE 86 : 1 CONTINUE 87 : ELSEIF (IESREF.EQ.2) THEN 88 : DO 3 IELEM=1,NDELEM 89 : DO 32 IPOGO=1,NDPOGO 90 : DET=(JMAJAC(1,1,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM)) 91 : $ - (JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(1,2,IPOGO,IELEM)) 92 : IF (DET.NE.ZERO) THEN 93 : IF (IPOGO.EQ.1) THEN 94 : LSIGN=(DET.GT.ZERO) 95 : ELSE 96 : IF (LSIGN.NEQV.(DET.GT.ZERO)) GOTO 9997 97 : ENDIF 98 : INVDET=UN/DET 99 : JDTJAC(IPOGO,IELEM)=DET 100 : JMIJAC(1,1,IPOGO,IELEM)= 101 : $ (+JMAJAC(2,2,IPOGO,IELEM))*INVDET 102 : JMIJAC(2,1,IPOGO,IELEM)= 103 : $ (-JMAJAC(2,1,IPOGO,IELEM))*INVDET 104 : JMIJAC(1,2,IPOGO,IELEM)= 105 : $ (-JMAJAC(1,2,IPOGO,IELEM))*INVDET 106 : JMIJAC(2,2,IPOGO,IELEM)= 107 : $ (+JMAJAC(1,1,IPOGO,IELEM))*INVDET 108 : ELSE 109 : GOTO 9998 110 : ENDIF 111 : 32 CONTINUE 112 : 3 CONTINUE 113 : ELSEIF (IESREF.EQ.3) THEN 114 : DO 5 IELEM=1,NDELEM 115 : DO 52 IPOGO=1,NDPOGO 116 : JMIJAC(1,1,IPOGO,IELEM)= 117 : $ (JMAJAC(2,2,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM)) 118 : $ -(JMAJAC(3,2,IPOGO,IELEM)*JMAJAC(2,3,IPOGO,IELEM)) 119 : JMIJAC(2,1,IPOGO,IELEM)= 120 : $ (JMAJAC(3,1,IPOGO,IELEM)*JMAJAC(2,3,IPOGO,IELEM)) 121 : $ -(JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM)) 122 : JMIJAC(3,1,IPOGO,IELEM)= 123 : $ (JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(3,2,IPOGO,IELEM)) 124 : $ -(JMAJAC(3,1,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM)) 125 : JMIJAC(1,2,IPOGO,IELEM)= 126 : $ (JMAJAC(1,3,IPOGO,IELEM)*JMAJAC(3,2,IPOGO,IELEM)) 127 : $ -(JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM)) 128 : JMIJAC(2,2,IPOGO,IELEM)= 129 : $ (JMAJAC(1,1,IPOGO,IELEM)*JMAJAC(3,3,IPOGO,IELEM)) 130 : $ -(JMAJAC(1,3,IPOGO,IELEM)*JMAJAC(3,1,IPOGO,IELEM)) 131 : JMIJAC(3,2,IPOGO,IELEM)= 132 : $ (JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(3,1,IPOGO,IELEM)) 133 : $ -(JMAJAC(3,2,IPOGO,IELEM)*JMAJAC(1,1,IPOGO,IELEM)) 134 : JMIJAC(1,3,IPOGO,IELEM)= 135 : $ (JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(2,3,IPOGO,IELEM)) 136 : $ -(JMAJAC(1,3,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM)) 137 : JMIJAC(2,3,IPOGO,IELEM)= 138 : $ (JMAJAC(2,1,IPOGO,IELEM)*JMAJAC(1,3,IPOGO,IELEM)) 139 : $ -(JMAJAC(2,3,IPOGO,IELEM)*JMAJAC(1,1,IPOGO,IELEM)) 140 : JMIJAC(3,3,IPOGO,IELEM)= 141 : $ (JMAJAC(1,1,IPOGO,IELEM)*JMAJAC(2,2,IPOGO,IELEM)) 142 : $ -(JMAJAC(1,2,IPOGO,IELEM)*JMAJAC(2,1,IPOGO,IELEM)) 143 : DET=(JMAJAC(1,1,IPOGO,IELEM)*JMIJAC(1,1,IPOGO,IELEM)) 144 : $ +(JMAJAC(1,2,IPOGO,IELEM)*JMIJAC(2,1,IPOGO,IELEM)) 145 : $ +(JMAJAC(1,3,IPOGO,IELEM)*JMIJAC(3,1,IPOGO,IELEM)) 146 : * IF (DET.NE.ZERO) THEN 147 : IF (ABS(DET).GT.(SQRT(XPETIT))) THEN 148 : IF (IPOGO.EQ.1) THEN 149 : LSIGN=(DET.GT.ZERO) 150 : ELSE 151 : IF (LSIGN.NEQV.(DET.GT.ZERO)) GOTO 9997 152 : ENDIF 153 : INVDET=UN/DET 154 : JDTJAC(IPOGO,IELEM)=DET 155 : DO 522 JREFER=1,IESREF 156 : DO 5222 IREFER=1,IESREF 157 : JMIJAC(IREFER,JREFER,IPOGO,IELEM)= 158 : $ JMIJAC(IREFER,JREFER,IPOGO,IELEM) 159 : $ * INVDET 160 : 5222 CONTINUE 161 : 522 CONTINUE 162 : ELSE 163 : GOTO 9998 164 : ENDIF 165 : 52 CONTINUE 166 : 5 CONTINUE 167 : ELSE 168 : WRITE(IOIMP,*) 'Je ne sais pas calculer l''inverse d''une' 169 : WRITE(IOIMP,*) 'matrice jacobienne de dimension ' 170 : WRITE(IOIMP,*) 'IESREF=',IESREF 171 : ENDIF 172 : * 173 : * Normal termination 174 : * 175 : IRET=0 176 : RETURN 177 : * 178 : * Format handling 179 : * 180 : * 181 : * Error handling 182 : * 183 : 9997 CONTINUE 184 : IF (LERJ) THEN 185 : IRET=666 185 : CALL ECRENT(IRET) 186 : RETURN 187 : ENDIF 188 : WRITE(IOIMP,*) 'Le déterminant de la matrice jacobienne' 189 : WRITE(IOIMP,*) 'change de signe sur l''élément.' 190 : WRITE(IOIMP,*) 'IELEM=',IELEM,' IPOGO=',IPOGO 191 : GOTO 9999 192 : 9998 CONTINUE 193 : IF (LERJ) THEN 194 : IRET=667 194 : CALL ECRENT(IRET) 195 : RETURN 196 : ENDIF 197 : WRITE(IOIMP,*) 'Déterminant de la matrice jacobienne nul' 198 : WRITE(IOIMP,*) 'IELEM=',IELEM,' IPOGO=',IPOGO 199 : GOTO 9999 200 : 9999 CONTINUE 201 : IRET=1 202 : WRITE(IOIMP,*) 'An error was detected in subroutine geoli2' 203 : RETURN 204 : * 205 : * End of subroutine GEOLI2 206 : * 207 : END 208 : 209 :
© Cast3M 2003 - All rights reserved.
Disclaimer