1 : C GEOLIK SOURCE GOUNAND 09/09/03 21:15:15 6481 2 : SUBROUTINE GEOLIK(DFFPG,JCOOR,NBELEM, 3 : $ JMAJAC,JMIJAC,JDTJAC,LERJ,IMREG, 4 : $ IMPR,IRET) 5 : IMPLICIT REAL*8 (A-H,O-Z) 6 : IMPLICIT INTEGER (I-N) 7 : C*********************************************************************** 8 : C NOM : GEOLIK 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 : PRCHVA (impression d'un segment de type MCHEVA) 45 : C GEOLI1 (calcul de la matrice jacobienne (fortran 77 46 : C )) 47 : C GEOLI2 (calcul de l'inverse et du déterminant de la 48 : C matrice jacobienne (fortran 77)) 49 : C GEOLI3 (calcul de sqrt(det(transpose(J)*J)), J 50 : C étant la matrice jacobienne (fortran 77)) 51 : C APPELE PAR : NLIN 52 : C*********************************************************************** 53 : C ENTREES : * DFFPG (type MCHEVA) : valeurs des dérivées 54 : C premières des fonctions d'interpolation de la 55 : C transformation géométrique aux points de gauss 56 : C sur l'élément de référence. 57 : C Structure (cf.include SMCHAEL) : 58 : C (1, nb. ddl, 1, dim.esp.réf, nb. poi. gauss, 1) 59 : C * JCOOR (type MCHEVA) : valeurs des ddl de la 60 : C transformation géométrique sur le maillage 61 : C élémentaire courant. 62 : C Structure (cf.include SMCHAEL) : 63 : C (1, nb. ddl, 1, dim. esp. réel, 64 : C 1, nb. éléments) 65 : C * NBELEM (type entier) : nombre d'éléments du 66 : C maillage élémentaire courant. 67 : C ENTREES/SORTIES : - 68 : C SORTIES : * JMAJAC (type MCHEVA) : valeurs de la matrice 69 : C jacobienne aux points de Gauss sur le maillage 70 : C élémentaire courant. 71 : C Structure (cf.include SMCHAEL) : 72 : C (1, 1, dim. esp. réel, dim. esp. référence, 73 : C nb. poi. gauss, nb. éléments) 74 : C * JMIJAC (type MCHEVA) : valeurs de l'inverse de 75 : C la matrice jacobienne aux points de Gauss sur 76 : C le maillage élémentaire courant. 77 : C Structure (cf.include SMCHAEL) : 78 : C (1, 1, dim. esp. référence, dim. esp. réel, 79 : C nb. poi. gauss, nb. éléments) 80 : C JMIJAC n'existe que si dim.esp.réf=dim.esp.réel 81 : C * JDTJAC (type MCHEVA) : valeurs du déterminant 82 : C de la matrice jacobienne aux points de Gauss 83 : C sur le maillage élémentaire courant. 84 : C Structure (cf.include SMCHAEL) : 85 : C (1, 1, 1, 1, nb. poi. gauss, nb. éléments) 86 : C Si la matrice jacobienne n'est pas carrée, on 87 : C a calculé sqrt(det(trans(J).J)) 88 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé 89 : C*********************************************************************** 90 : C VERSION : v1, 09/08/99, version initiale 91 : C HISTORIQUE : v1, 09/08/99, création 92 : C HISTORIQUE : 93 : C HISTORIQUE : 94 : C*********************************************************************** 95 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 96 : C en cas de modification de ce sous-programme afin de faciliter 97 : C la maintenance ! 98 : C*********************************************************************** 99 : -INC CCOPTIO 100 : CBEGININCLUDE SMCHAEL 101 : SEGMENT MCHAEL 102 : POINTEUR IMACHE(N1).MELEME 103 : POINTEUR ICHEVA(N1).MCHEVA 104 : ENDSEGMENT 105 : SEGMENT MCHEVA 106 : REAL*8 VELCHE(NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM) 107 : ENDSEGMENT 108 : SEGMENT LCHEVA 109 : POINTEUR LISCHE(NBCHE).MCHEVA 110 : ENDSEGMENT 111 : CENDINCLUDE SMCHAEL 112 : INTEGER NBLIG,NBCOL,N2LIG,N2COL,NBPOI,NBELM 113 : POINTEUR JCOOR.MCHEVA 114 : POINTEUR JMAJAC.MCHEVA 115 : POINTEUR JMIJAC.MCHEVA 116 : POINTEUR JDTJAC.MCHEVA 117 : POINTEUR JJTJ.MCHEVA 118 : POINTEUR JJTJM1.MCHEVA 119 : * Valeur des fns d'interpolation de l'espace géométrique aux pts de Gauss 120 : POINTEUR DFFPG.MCHEVA 121 : * 122 : INTEGER NBELEM 123 : INTEGER IMPR,IRET 124 : * 125 : LOGICAL LCARRE,LERJ 126 : INTEGER IESREF,IESREL,NDNOEU,NDCOL,NDPOGO,NDELEM 127 : * 128 : * Executable statements 129 : * 130 : IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans geolik' 131 : SEGACT DFFPG 132 : NDNOEU=DFFPG.VELCHE(/2) 133 : IESREF=DFFPG.VELCHE(/4) 134 : NDPOGO=DFFPG.VELCHE(/5) 135 : SEGACT JCOOR 136 : NDCOL=JCOOR.VELCHE(/2) 137 : IESREL=JCOOR.VELCHE(/4) 138 : NDELEM=JCOOR.VELCHE(/6) 139 : IF (NDCOL.NE.NDNOEU) THEN 140 : WRITE(IOIMP,*) 'Incompatibilité fns d''interpolation-géométrie' 141 : WRITE(IOIMP,*) 'NDNOEU=',NDNOEU,' NDCOL=',NDCOL 142 : GOTO 9999 143 : ENDIF 144 : * On pourrait envisager un cas où tous les éléments seraient 145 : * identiques (ex. carré) NDELEM=1 146 : IF (NDELEM.NE.NBELEM) THEN 147 : WRITE(IOIMP,*) 'Incompatibilité coordonnées-géométrie' 148 : WRITE(IOIMP,*) 'NDELEM=',NDELEM,' NBELEM=',NBELEM 149 : GOTO 9999 150 : ENDIF 151 : IF (IMREG.EQ.1) THEN 152 : NDELEM=1 153 : NDPOGO=1 154 : ENDIF 155 : IF (IESREF.LT.IESREL) THEN 156 : LCARRE=.FALSE. 157 : ELSEIF (IESREF.EQ.IESREL) THEN 158 : LCARRE=.TRUE. 159 : ELSE 160 : WRITE(IOIMP,*) 'Dim. esp. ref. > dim. esp. reel ???' 161 : WRITE(IOIMP,*) 'IESREF=',IESREF,' IESREL=',IESREL 162 : GOTO 9999 163 : ENDIF 164 : * 165 : * Initialisations... 166 : * 167 : NBLIG=1 168 : NBCOL=1 169 : N2LIG=IESREL 170 : N2COL=IESREF 171 : NBPOI=NDPOGO 172 : NBELM=NDELEM 173 : SEGINI JMAJAC 174 : * 175 : * On effectue le calcul de la matrice jacobienne aux points de Gauss 176 : * 176 : CALL GEOLI1(IESREL,IESREF,NDNOEU,NDPOGO,NDELEM, 177 : $ DFFPG.VELCHE,JCOOR.VELCHE, 178 : $ JMAJAC.VELCHE, 179 : $ IMPR,IRET) 180 : IF (IRET.NE.0) GOTO 9999 181 : IF (IMPR.GT.4) THEN 182 : WRITE(IOIMP,*) 'On a créé', 183 : $ ' JMAJAC(élément , poi.gauss ,', 184 : $ ' dksi(i) , dx(j) ,1,1)' 184 : CALL PRCHVA(JMAJAC,IMPR,IRET) 185 : IF (IRET.NE.0) GOTO 9999 186 : ENDIF 187 : * 188 : * Si la matrice est carrée, on tente de calculer 189 : * son inverse et son déterminant... 190 : * 191 : IF (LCARRE) THEN 192 : NBLIG=1 193 : NBCOL=1 194 : N2LIG=IESREF 195 : N2COL=IESREL 196 : NBPOI=NDPOGO 197 : NBELM=NDELEM 198 : SEGINI JMIJAC 199 : NBLIG=1 200 : NBCOL=1 201 : N2LIG=1 202 : N2COL=1 203 : NBPOI=NDPOGO 204 : NBELM=NDELEM 205 : SEGINI JDTJAC 205 : CALL GEOLI2(IESREF,NDPOGO,NDELEM,JMAJAC.VELCHE, 206 : $ JMIJAC.VELCHE,JDTJAC.VELCHE,LERJ, 207 : $ IMPR,IRET) 208 : IF (IRET.NE.0) THEN 209 : IF (LERJ) GOTO 9666 210 : GOTO 9999 211 : ENDIF 212 : SEGDES JDTJAC 213 : SEGDES JMIJAC 214 : IF (IMPR.GT.4) THEN 215 : WRITE(IOIMP,*) 'On a créé', 216 : $ ' JMIJAC(élément , poi.gauss ,', 217 : $ ' dx(i) , dksi(j) ,1,1)' 217 : CALL PRCHVA(JMIJAC,IMPR,IRET) 218 : IF (IRET.NE.0) GOTO 9999 219 : ENDIF 220 : ELSE 221 : * 222 : * Sinon on calcule : SQRT ( det (transpose(A) * A)) 223 : * et la pseudo-inverse (JtJ)^-1 Jt 224 : * 225 : NBLIG=1 226 : NBCOL=1 227 : N2LIG=IESREF 228 : N2COL=IESREL 229 : NBPOI=NDPOGO 230 : NBELM=NDELEM 231 : SEGINI JMIJAC 232 : * 233 : NBLIG=1 234 : NBCOL=1 235 : N2LIG=1 236 : N2COL=1 237 : NBPOI=NDPOGO 238 : NBELM=NDELEM 239 : SEGINI JDTJAC 240 : * Objets temporaire 241 : NBLIG=1 242 : NBCOL=1 243 : N2LIG=IESREF 244 : N2COL=IESREF 245 : NBPOI=NDPOGO 246 : NBELM=NDELEM 247 : SEGINI JJTJ 248 : SEGINI JJTJM1 248 : CALL GEOLI4(IESREL,IESREF,NDPOGO,NDELEM,JMAJAC.VELCHE, 249 : $ JJTJ.VELCHE,JJTJM1.VELCHE, 250 : $ JMIJAC.VELCHE,JDTJAC.VELCHE, 251 : $ IMPR,IRET) 252 : IF (IRET.NE.0) GOTO 9999 253 : SEGSUP JJTJM1 254 : SEGSUP JJTJ 255 : SEGDES JMIJAC 256 : SEGDES JDTJAC 257 : ENDIF 258 : SEGDES JMAJAC 259 : SEGDES JCOOR 260 : SEGDES DFFPG 261 : IF (IMPR.GT.3) THEN 262 : * IF (IMREG.GT.0) THEN 263 : WRITE(IOIMP,*) 'On a créé', 264 : $ ' JDTJAC(élément , poi.gauss ,', 265 : $ '1,1,1,1)' 265 : CALL PRCHVA(JDTJAC,IMPR,IRET) 265 : * CALL PRCHVA(JDTJAC,6,IRET) 266 : IF (IRET.NE.0) GOTO 9999 267 : ENDIF 268 : * ENDIF 269 : * 270 : * Normal termination 271 : * 272 : IRET=0 273 : RETURN 274 : * 275 : * Format handling 276 : * 277 : * 278 : * Error handling 279 : * 280 : 9666 CONTINUE 281 : IRET=666 282 : RETURN 283 : 9999 CONTINUE 284 : IRET=1 285 : WRITE(IOIMP,*) 'An error was detected in subroutine geolik' 286 : RETURN 287 : * 288 : * End of subroutine GEOLIK 289 : * 290 : END 291 : 292 : 293 : 294 :
© Cast3M 2003 - All rights reserved.
Disclaimer