1 : C GEOLIF SOURCE GOUNAND 06/08/04 21:15:57 5520 2 : SUBROUTINE GEOLIF(DFFPG,JCOOR,SSFACT,NBELEF, 3 : $ JMAJAC,JMIJAC,JDTJAC,LERJ, 4 : $ IMPR,IRET) 5 : IMPLICIT REAL*8 (A-H,O-Z) 6 : IMPLICIT INTEGER (I-N) 7 : C*********************************************************************** 8 : C NOM : GEOLIF 9 : C PROJET : Noyau linéaire NLIN 10 : C DESCRIPTION : Calcul de la matrice jacobienne d'une transformation 11 : C géométrique aux points de Gauss d'un élément de 12 : C référence pour chaque élément réel. 13 : C On a : 14 : C Fonction f : R^m -> R^n 15 : C a=(a1...am) |-> f(a)=(f1(a)...fn(a) 16 : C => matjac(i,j)(a)=Dj fi (a) = dfi / dxj (a) 17 : C 18 : C Par exemple, pour une surface en 3D, la matrice jacobienne 19 : C s'exprime comme : 20 : C 21 : C [] 22 : C [J] = [ ] . [ { dNg/d(ksi) } { dNg/d(eta) } ] 23 : C [ ] 24 : C 25 : C (3x2) (3 x Nnoeuds) (Nnoeuds x 2) 26 : C 27 : C Ici, le nb de ddl est égal aux nbs de noeuds car l'interpolation 28 : C pour la géométrie est de type nodale. 29 : C 30 : C Si la matrice jacobienne J est carrée, on calcule 31 : C également l'inverse de la matrice jacobienne j ainsi que 32 : C son déterminant : det J (appelé Jacobien). 33 : C Si la matrice jacobienne n'est pas carrée, le Jacobien 34 : C est calculé comme : 35 : C sqrt(det(transpose(J)*J)) 36 : C Le but de ce sous-programme est d'effectuer toutes les 37 : C opérations de calcul dépendant uniquement de la 38 : C géométrie (ie de la variable : coordonnées des points). 39 : C 40 : C LANGAGE : ESOPE 41 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF) 42 : C mél : gounand@semt2.smts.cea.fr 43 : C*********************************************************************** 44 : C APPELES : RSET8 (copie de tableaux de CHARACTER*8) 45 : C PRCHVA (impression d'un segment de type MCHEVA) 46 : C GEOLI1 (calcul de la matrice jacobienne (fortran 77 47 : C )) 48 : C GEOLI2 (calcul de l'inverse et du déterminant de la 49 : C matrice jacobienne (fortran 77)) 50 : C GEOLI3 (calcul de sqrt(det(transpose(J)*J)), J 51 : C étant la matrice jacobienne (fortran 77)) 52 : C APPELE PAR : NLIN 53 : C*********************************************************************** 54 : C ENTREES : * DFFPG (type MCHEVA) : valeurs des dérivées 55 : C premières des fonctions d'interpolation de la 56 : C transformation géométrique aux points de gauss 57 : C sur l'élément de référence. 58 : C Structure (cf.include SMCHAEL) : 59 : C (1, nb. ddl, 1, dim.esp.réf, nb. poi. gauss, 1) 60 : C * JCOOR (type MCHEVA) : valeurs des ddl de la 61 : C transformation géométrique sur le maillage 62 : C élémentaire courant. 63 : C Structure (cf.include SMCHAEL) : 64 : C (1, nb. ddl, 1, dim. esp. réel, 65 : C 1, nb. éléments) 66 : C * NBELEV (type entier) : nombre d'éléments du 67 : C maillage élémentaire courant. 68 : C ENTREES/SORTIES : - 69 : C SORTIES : * JMAJAC (type MCHEVA) : valeurs de la matrice 70 : C jacobienne aux points de Gauss sur le maillage 71 : C élémentaire courant. 72 : C Structure (cf.include SMCHAEL) : 73 : C (1, 1, dim. esp. réel, dim. esp. référence, 74 : C nb. poi. gauss, nb. éléments) 75 : C * JMIJAC (type MCHEVA) : valeurs de l'inverse de 76 : C la matrice jacobienne aux points de Gauss sur 77 : C le maillage élémentaire courant. 78 : C Structure (cf.include SMCHAEL) : 79 : C (1, 1, dim. esp. référence, dim. esp. réel, 80 : C nb. poi. gauss, nb. éléments) 81 : C JMIJAC n'existe que si dim.esp.réf=dim.esp.réel 82 : C * JDTJAC (type MCHEVA) : valeurs du déterminant 83 : C de la matrice jacobienne aux points de Gauss 84 : C sur le maillage élémentaire courant. 85 : C Structure (cf.include SMCHAEL) : 86 : C (1, 1, 1, 1, nb. poi. gauss, nb. éléments) 87 : C Si la matrice jacobienne n'est pas carrée, on 88 : C a calculé sqrt(det(trans(J).J)) 89 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé 90 : C*********************************************************************** 91 : C VERSION : v2, 03/10/03, refonte complète (modif SMPOUET) 92 : C VERSION : v1, 09/08/99, version initiale 93 : C HISTORIQUE : v1, 09/08/99, création 94 : C HISTORIQUE : 95 : C HISTORIQUE : 96 : C*********************************************************************** 97 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 98 : C en cas de modification de ce sous-programme afin de faciliter 99 : C la maintenance ! 100 : C*********************************************************************** 101 : -INC CCOPTIO 102 : CBEGININCLUDE SFACTIV 103 : SEGMENT FACTIV 104 : POINTEUR IFACTI(NBSOUV).SFACTI 105 : ENDSEGMENT 106 : SEGMENT SFACTI 107 : POINTEUR ISFACT(NBSOFV).SSFACT 108 : ENDSEGMENT 109 : SEGMENT SSFACT 110 : LOGICAL LFACTI(NBELFV,NBELEV) 111 : ENDSEGMENT 112 : CENDINCLUDE SFACTIV 113 : CBEGININCLUDE SMCHAEL 114 : SEGMENT MCHAEL 115 : POINTEUR IMACHE(N1).MELEME 116 : POINTEUR ICHEVA(N1).MCHEVA 117 : ENDSEGMENT 118 : SEGMENT MCHEVA 119 : REAL*8 VELCHE(NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM) 120 : ENDSEGMENT 121 : SEGMENT LCHEVA 122 : POINTEUR LISCHE(NBCHE).MCHEVA 123 : ENDSEGMENT 124 : CENDINCLUDE SMCHAEL 125 : INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM 126 : POINTEUR JCOOR.MCHEVA 127 : POINTEUR JMAJAC.MCHEVA 128 : POINTEUR JMIJAC.MCHEVA 129 : POINTEUR JDTJAC.MCHEVA 130 : * Valeur des fns d'interpolation de l'espace géométrique aux pts de Gauss 131 : POINTEUR DFFPG.MCHEVA 132 : * 133 : INTEGER NBELEV,NBELEF,NBELFV 134 : INTEGER IMPR,IRET 135 : * 136 : LOGICAL LCARRE,LERJ 137 : INTEGER IESREF,IESREL,NBPOGO 138 : * 139 : * Executable statements 140 : * 141 : IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans geolif' 142 : SEGACT SSFACT 143 : NBELFV=SSFACT.LFACTI(/1) 144 : NBELEV=SSFACT.LFACTI(/2) 145 : SEGACT JCOOR 146 : NBLIG=JCOOR.VELCHE(/1) 147 : NBCOL=JCOOR.VELCHE(/2) 148 : N2LIG=JCOOR.VELCHE(/3) 149 : N2COL=JCOOR.VELCHE(/4) 150 : NBPOI=JCOOR.VELCHE(/5) 151 : NBELM=JCOOR.VELCHE(/6) 152 : * On pourrait envisager un cas où tous les éléments seraient 153 : * identiques (ex. carré) NBELM=1 154 : IF (NBLIG.NE.1.OR.N2LIG.NE.1 155 : $ .OR.NBPOI.NE.1.OR.NBELM.NE.NBELEV) THEN 156 : WRITE(IOIMP,*) 'Erreur dims JCOOR' 157 : GOTO 9999 158 : ENDIF 159 : NDDL=NBCOL 160 : IESREL=N2COL 161 : SEGACT DFFPG 162 : NBLIG=DFFPG.VELCHE(/1) 163 : NBCOL=DFFPG.VELCHE(/2) 164 : N2LIG=DFFPG.VELCHE(/3) 165 : N2COL=DFFPG.VELCHE(/4) 166 : NBPOI=DFFPG.VELCHE(/5) 167 : NBELM=DFFPG.VELCHE(/6) 168 : IF (NBLIG.NE.1.OR.NBCOL.NE.NDDL.OR.N2LIG.NE.1 169 : $ .OR.(NBELM.NE.NBELFV.AND.NBELM.NE.1)) THEN 170 : WRITE(IOIMP,*) 'Erreur dims DFFPG' 171 : GOTO 9999 172 : ENDIF 173 : IESREF=N2COL 174 : NBPOGO=NBPOI 175 : NLFVDF=NBELM 176 : IF (IESREF.LT.IESREL) THEN 177 : LCARRE=.FALSE. 178 : ELSEIF (IESREF.EQ.IESREL) THEN 179 : LCARRE=.TRUE. 180 : ELSE 181 : WRITE(IOIMP,*) 'Dim. esp. ref. > dim. esp. reel ???' 182 : WRITE(IOIMP,*) 'IESREF=',IESREF,' IESREL=',IESREL 183 : GOTO 9999 184 : ENDIF 185 : * 186 : * Initialisations... 187 : * 188 : NBLIG=1 189 : NBCOL=1 190 : N2LIG=IESREL 191 : N2COL=IESREF 192 : NBPOI=NBPOGO 193 : NBELM=NBELEF 194 : SEGINI JMAJAC 195 : * SEGPRT,DFFPG 196 : * 197 : * On effectue le calcul de la matrice jacobienne aux points de Gauss 198 : * 198 : CALL GEOLF1(IESREL,IESREF,NDDL,NBPOGO,NBELEV,NBELFV,NBELEF, 199 : $ NLFVDF, 200 : $ DFFPG.VELCHE,JCOOR.VELCHE,SSFACT.LFACTI, 201 : $ JMAJAC.VELCHE, 202 : $ IMPR,IRET) 203 : IF (IRET.NE.0) GOTO 9999 204 : * SEGPRT,JMAJAC 205 : * 206 : * Si la matrice est carrée, on tente de calculer 207 : * son inverse et son déterminant... 208 : * 209 : IF (LCARRE) THEN 210 : NBLIG=1 211 : NBCOL=1 212 : N2LIG=IESREF 213 : N2COL=IESREL 214 : NBPOI=NBPOGO 215 : NBELM=NBELEF 216 : SEGINI JMIJAC 217 : NBLIG=1 218 : NBCOL=1 219 : N2LIG=1 220 : N2COL=1 221 : NBPOI=NBPOGO 222 : NBELM=NBELEF 223 : SEGINI JDTJAC 223 : CALL GEOLI2(IESREF,NBPOGO,NBELEF,JMAJAC.VELCHE, 224 : $ JMIJAC.VELCHE,JDTJAC.VELCHE,LERJ, 225 : $ IMPR,IRET) 226 : * SEGPRT,JDTJAC 227 : IF (IRET.NE.0) THEN 228 : IF (LERJ) GOTO 9666 229 : GOTO 9999 230 : ENDIF 231 : SEGDES JDTJAC 232 : SEGDES JMIJAC 233 : ELSE 234 : * 235 : * Sinon on calcule : SQRT ( det (transpose(A) * A)) 236 : * 237 : * Pas de matrice inverse 238 : JMIJAC=0 239 : * 240 : NBLIG=1 241 : NBCOL=1 242 : N2LIG=1 243 : N2COL=1 244 : NBPOI=NBPOGO 245 : NBELM=NBELEF 246 : SEGINI JDTJAC 246 : CALL GEOLI3(IESREL,IESREF,NBPOGO,NBELEF,JMAJAC.VELCHE, 247 : $ JDTJAC.VELCHE, 248 : $ IMPR,IRET) 249 : IF (IRET.NE.0) GOTO 9999 250 : SEGDES JDTJAC 251 : ENDIF 252 : SEGDES JMAJAC 253 : SEGDES SSFACT 254 : SEGDES JCOOR 255 : SEGDES DFFPG 256 : * 257 : * Normal termination 258 : * 259 : IRET=0 260 : RETURN 261 : * 262 : * Format handling 263 : * 264 : * 265 : * Error handling 266 : * 267 : 9666 CONTINUE 268 : IRET=666 269 : RETURN 270 : 9999 CONTINUE 271 : IRET=1 272 : WRITE(IOIMP,*) 'An error was detected in subroutine geolif' 273 : RETURN 274 : * 275 : * End of subroutine GEOLIF 276 : * 277 : END 278 : 279 :
© Cast3M 2003 - All rights reserved.
Disclaimer