1 : C GEOLIN SOURCE GOUNAND 06/11/24 21:15:15 5574 2 : SUBROUTINE GEOLIN(DFFPG,JCOOR,NBELEM, 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 : GEOLIN 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 geolin' 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 (IESREF.LT.IESREL) THEN 152 : LCARRE=.FALSE. 153 : ELSEIF (IESREF.EQ.IESREL) THEN 154 : LCARRE=.TRUE. 155 : ELSE 156 : WRITE(IOIMP,*) 'Dim. esp. ref. > dim. esp. reel ???' 157 : WRITE(IOIMP,*) 'IESREF=',IESREF,' IESREL=',IESREL 158 : GOTO 9999 159 : ENDIF 160 : * 161 : * Initialisations... 162 : * 163 : NBLIG=1 164 : NBCOL=1 165 : N2LIG=IESREL 166 : N2COL=IESREF 167 : NBPOI=NDPOGO 168 : NBELM=NDELEM 169 : SEGINI JMAJAC 170 : * 171 : * On effectue le calcul de la matrice jacobienne aux points de Gauss 172 : * 172 : CALL GEOLI1(IESREL,IESREF,NDNOEU,NDPOGO,NDELEM, 173 : $ DFFPG.VELCHE,JCOOR.VELCHE, 174 : $ JMAJAC.VELCHE, 175 : $ IMPR,IRET) 176 : IF (IRET.NE.0) GOTO 9999 177 : IF (IMPR.GT.4) THEN 178 : WRITE(IOIMP,*) 'On a créé', 179 : $ ' JMAJAC(élément , poi.gauss ,', 180 : $ ' dksi(i) , dx(j) ,1,1)' 180 : CALL PRCHVA(JMAJAC,IMPR,IRET) 181 : IF (IRET.NE.0) GOTO 9999 182 : ENDIF 183 : * 184 : * Si la matrice est carrée, on tente de calculer 185 : * son inverse et son déterminant... 186 : * 187 : IF (LCARRE) THEN 188 : NBLIG=1 189 : NBCOL=1 190 : N2LIG=IESREF 191 : N2COL=IESREL 192 : NBPOI=NDPOGO 193 : NBELM=NDELEM 194 : SEGINI JMIJAC 195 : NBLIG=1 196 : NBCOL=1 197 : N2LIG=1 198 : N2COL=1 199 : NBPOI=NDPOGO 200 : NBELM=NDELEM 201 : SEGINI JDTJAC 201 : CALL GEOLI2(IESREF,NDPOGO,NDELEM,JMAJAC.VELCHE, 202 : $ JMIJAC.VELCHE,JDTJAC.VELCHE,LERJ, 203 : $ IMPR,IRET) 204 : IF (IRET.NE.0) THEN 205 : IF (LERJ) GOTO 9666 206 : GOTO 9999 207 : ENDIF 208 : SEGDES JDTJAC 209 : SEGDES JMIJAC 210 : IF (IMPR.GT.4) THEN 211 : WRITE(IOIMP,*) 'On a créé', 212 : $ ' JMIJAC(élément , poi.gauss ,', 213 : $ ' dx(i) , dksi(j) ,1,1)' 213 : CALL PRCHVA(JMIJAC,IMPR,IRET) 214 : IF (IRET.NE.0) GOTO 9999 215 : ENDIF 216 : ELSE 217 : * 218 : * Sinon on calcule : SQRT ( det (transpose(A) * A)) 219 : * et la pseudo-inverse (JtJ)^-1 Jt 220 : * 221 : NBLIG=1 222 : NBCOL=1 223 : N2LIG=IESREF 224 : N2COL=IESREL 225 : NBPOI=NDPOGO 226 : NBELM=NDELEM 227 : SEGINI JMIJAC 228 : * 229 : NBLIG=1 230 : NBCOL=1 231 : N2LIG=1 232 : N2COL=1 233 : NBPOI=NDPOGO 234 : NBELM=NDELEM 235 : SEGINI JDTJAC 236 : * Objets temporaire 237 : NBLIG=1 238 : NBCOL=1 239 : N2LIG=IESREF 240 : N2COL=IESREF 241 : NBPOI=NDPOGO 242 : NBELM=NDELEM 243 : SEGINI JJTJ 244 : SEGINI JJTJM1 244 : CALL GEOLI4(IESREL,IESREF,NDPOGO,NDELEM,JMAJAC.VELCHE, 245 : $ JJTJ.VELCHE,JJTJM1.VELCHE, 246 : $ JMIJAC.VELCHE,JDTJAC.VELCHE, 247 : $ IMPR,IRET) 248 : IF (IRET.NE.0) GOTO 9999 249 : SEGSUP JJTJM1 250 : SEGSUP JJTJ 251 : SEGDES JMIJAC 252 : SEGDES JDTJAC 253 : ENDIF 254 : SEGDES JMAJAC 255 : SEGDES JCOOR 256 : SEGDES DFFPG 257 : IF (IMPR.GT.3) THEN 258 : WRITE(IOIMP,*) 'On a créé', 259 : $ ' JDTJAC(élément , poi.gauss ,', 260 : $ '1,1,1,1)' 260 : CALL PRCHVA(JDTJAC,IMPR,IRET) 261 : IF (IRET.NE.0) GOTO 9999 262 : ENDIF 263 : * 264 : * Normal termination 265 : * 266 : IRET=0 267 : RETURN 268 : * 269 : * Format handling 270 : * 271 : * 272 : * Error handling 273 : * 274 : 9666 CONTINUE 275 : IRET=666 276 : RETURN 277 : 9999 CONTINUE 278 : IRET=1 279 : WRITE(IOIMP,*) 'An error was detected in subroutine geolin' 280 : RETURN 281 : * 282 : * End of subroutine GEOLIN 283 : * 284 : END 285 : 286 : 287 :
© Cast3M 2003 - All rights reserved.
Disclaimer