1 : C MELIM SOURCE MAGN 13/02/06 21:15:01 7695 2 : SUBROUTINE MELIM(MATRIK,KCLIM, 3 : $ INCX,KS2B, 4 : $ KMORS,KISA, 5 : $ IMPR,IRET) 6 : IMPLICIT REAL*8 (A-H,O-Z) 7 : IMPLICIT INTEGER (I-N) 8 : C*********************************************************************** 9 : C NOM : MELIM 10 : C DESCRIPTION : 11 : C Ce SP s'occupe de prendre en compte les conditions 12 : C aux limites de type Dirichlet pour un système linéaire 13 : C du type Ax=b (avec A utilisant le stockage Morse. 14 : C On procède en trois étapes : 15 : C 1 - on surcharge le second membre (b) 16 : C et l'inconnue (x) avec le chpoint des conditions 17 : C aux limites 18 : C 2 - on s'occupe de la matrice Morse (A) 19 : C Aux degrés de libertés où il y a des 20 : C conditions aux limites : 1 sur la diagonale 21 : C 0 sur la ligne 22 : C 0 sur la colonne en reportant les 23 : C produits dans le second membre 24 : C (sinon perte de symétrie 25 : C éventuelle de A) 26 : C 3 - On nettoie la matrice morse A des 0 qu'on vient de lui 27 : C rajouter (sous-programme clmors). 28 : C 29 : C 2 et 3 peuvent sans doute etre regroupées au prix d'une 30 : C perte de lisibilité 31 : C 32 : C 33 : C LANGAGE : ESOPE 34 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF) 35 : C mél : gounand@semt2.smts.cea.fr 36 : C*********************************************************************** 37 : C APPELES : KRIPAD, CLMORS 38 : C*********************************************************************** 39 : C ENTREES : MATRIK, KCLIM, IMPR 40 : C ENTREES/SORTIES : KS2B, INCX 41 : C SORTIES : KMORS,KISA,IRET 42 : C CODE RETOUR (IRET) : 0 si ok 43 : C <0 si problème 44 : C MATRIK : pointeur sur segment MATRIK de l'include SMMATRIK 45 : C on pioche dedans les informations nécessaires 46 : C (différents pointeurs, nb. de ddl...) 47 : C KCLIM : pointeur sur segment MCHPOI de l'include SMCHPOI 48 : C chpoint contenant les conditions aux limites de 49 : C Dirichlet (valeur imposée). 50 : C IMPR : niveau d'impression 51 : C KS2B : pointeur sur segment IZA de l'include SMMATRIK 52 : C vecteur second membre 53 : C (ie b dans la résolution de Ax=b). 54 : C INCX : pointeur sur segment IZA de l'include SMMATRIK 55 : C vecteur des inconnues 56 : C (ie x dans la résolution de Ax=b). 57 : C KMORS : pointeurs sur segment PMORS et IZA de l'include 58 : C KISA SMMATRIK : profil et valeurs de la matrice morse 59 : C (ie A dans la résolution de Ax=b). 60 : C*********************************************************************** 61 : C VERSION : v1, 20/12/99 62 : C HISTORIQUE : v1, 01/04/98, création 63 : C HISTORIQUE : 20/12/99, 64 : C Nouvelle signification de NUAN (renumérotation des ddl et non pas 65 : C des points) 66 : C HISTORIQUE : 05/01/00 : non-duplication de la matrice assemblée 67 : C HISTORIQUE : 09/04/04 : idmatp idmatd, 68 : C la cl dirichlet n'est plus forcément sur la 69 : C diagonale 70 : C HISTORIQUE : 71 : C*********************************************************************** 72 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 73 : C en cas de modification de ce sous-programme afin de faciliter 74 : C la maintenance ! 75 : C*********************************************************************** 76 : -INC CCOPTIO 77 : -INC SMELEME 78 : POINTEUR ISPGT.MELEME 79 : -INC SMLENTI 80 : POINTEUR DDLPD.MLENTI 81 : POINTEUR DDLDP.MLENTI 82 : *-INC SMLLOGI 83 : * SEGMENT MLLOGI 84 : * LOGICAL LOGI(JG) 85 : * ENDSEGMENT 86 : * INTEGER JG 87 : * POINTEUR DDLCL.MLLOGI 88 : -INC SMCHPOI 89 : POINTEUR KCLIM.MCHPOI 90 : -INC SMMATRIK 91 : POINTEUR KMORS.PMORS 92 : POINTEUR KISA.IZA 93 : POINTEUR AMORS.PMORS 94 : POINTEUR AISA.IZA 95 : POINTEUR KS2B.IZA 96 : POINTEUR INCX.IZA 97 : POINTEUR IDMATP.IDMAT 98 : POINTEUR IDMATD.IDMAT 99 : * 100 : INTEGER IMPR,IRET 101 : * 102 : * .. Variables locales 103 : * .. Parameters 104 : REAL*8 ZERO,ONE 105 : PARAMETER (ZERO=0.0D0, ONE=1.0D0) 106 : * .. 107 : CHARACTER*8 NOMINC 108 : LOGICAL FLINC 109 : LOGICAL FLDIA 110 : * 111 : INTEGER I1,IN,ISOUPO,INC,INBI,INBVA 112 : INTEGER N ,NSOUPO,NC ,NBI 113 : INTEGER ICOL,INTTT,JNTTT 114 : INTEGER NTTT 115 : REAL*8 TEMP 116 : LOGICAL LINCX 117 : C*** 118 : LINCX=(INCX.NE.0) 119 : C On récupère les segments utiles 120 : IF (IMPR.GT.5) THEN 121 : WRITE(IOIMP,*) 'melim.eso : CL de Dirichlet' 122 : ENDIF 123 : SEGACT MATRIK 124 : MINC =KMINC 125 : ISPGT=KISPGT 126 : NTTT =KNTTT 127 : IDMATP=KIDMAT(1) 128 : IDMATD=KIDMAT(2) 129 : AMORS=KIDMAT(4) 130 : AISA =KIDMAT(5) 131 : SEGDES MATRIK 132 : C S'il n'y a pas de CL de Dirichlet, on n'a rien à faire 133 : IF(KCLIM.EQ.0) THEN 134 : SEGINI,KMORS=AMORS 135 : SEGDES KMORS 136 : SEGINI,KISA=AISA 137 : SEGDES KISA 138 : ELSE 139 : C 1) On s'occupe du second membre et de l'inconnue 140 : C On parcourt le chpoint un peu comme dans ch2vec.eso 141 : C Les DDL ou on applique une condition de Dirichlet 142 : C ont un DDLCL(No DDL)=1 et 0 sinon 143 : IF (IMPR.GT.6) THEN 144 : WRITE(IOIMP,*) 145 : $ 'Dirichlet sur l''inconnue et le second membre' 146 : ENDIF 147 : JG=NTTT 148 : SEGINI DDLPD 149 : SEGINI DDLDP 150 : SEGACT KS2B*MOD 151 : IF (LINCX) SEGACT INCX*MOD 152 : MCHPOI=KCLIM 153 : SEGACT MINC 154 : NBI=LISINC(/2) 155 : C In KRIPAD : SEGACT ISPGT 156 : C SEGINI MLENTI 156 : CALL KRIPAD(ISPGT,MLENTI) 157 : SEGDES ISPGT 158 : SEGACT IDMATP 159 : SEGACT IDMATD 160 : SEGACT MCHPOI 161 : NSOUPO=IPCHP(/1) 162 : DO 1 ISOUPO=1,NSOUPO 163 : MSOUPO=IPCHP(ISOUPO) 164 : SEGACT MSOUPO 165 : NC=NOCOMP(/2) 166 : MELEME=IGEOC 167 : MPOVAL=IPOVAL 168 : SEGACT MELEME 169 : SEGACT MPOVAL 170 : C NPTI=VPOCHA(/1) 171 : DO 11 INC=1,NC 172 : NOMINC=NOCOMP(INC)//' ' 173 : FLINC=.FALSE. 174 : C Repeat..until 175 : INBI=1 176 : 111 CONTINUE 177 : IF (NOMINC.EQ.LISINC(INBI)) THEN 178 : FLINC=.TRUE. 179 : ELSEIF (INBI.LT.NBI) THEN 180 : INBI=INBI+1 181 : GOTO 111 182 : ENDIF 183 : IF (.NOT.FLINC) THEN 184 : * WRITE(IOIMP,*) 'melim.eso : composante ',NOMINC 185 : * $ ,'unknown' 186 : WRITE(IOIMP,*) 'CLIM : component name ',NOMINC 187 : $ ,' not in the matrix' 188 : GOTO 9999 189 : ELSE 190 : N=VPOCHA(/1) 191 : DO 112 IN=1,N 192 : I1=LECT(NUM(1,IN)) 193 : IF(I1.EQ.0)THEN 194 : * WRITE(IOIMP,*) 'melim.eso : le point ',IN 195 : * $ ,' du ch.' 196 : * WRITE(IOIMP,*) 'n''appartient pas au vec.' 197 : WRITE(IOIMP,*) 'CLIM : point number ',(num(1,IN)), 198 : $ ' unknown ',NOMINC,' not in the matrix' 199 : GOTO 9999 200 : ELSE 201 : IF (MPOS(I1,INBI).EQ.0) THEN 202 : WRITE(IOIMP,*) 'CLIM : point number ',(num(1,IN)), 203 : $ ' unknown ',NOMINC,' not in the matrix' 204 : GOTO 9999 205 : ELSE 206 : IDIAGP=IDMATP.NUAN(NPOS(I1)+MPOS(I1,INBI)-1) 207 : IDIAGD=IDMATD.NUAN(NPOS(I1)+MPOS(I1,INBI)-1) 208 : C* INBVA=NUAN(NPOS(I1)+MPOS(I1,INBI)-1) 209 : C INPT=NUAN(I1) 210 : C INBVA=NPOS(INPT)+MPOS(INPT,INBI)-1 211 : C KVEC.A(INBVA)=KVEC.A(INBVA)+VPOCHA(IN,INC) 212 : C Je préfère surcharger 213 : C 214 : * KS2B.A(INBVA)=VPOCHA(IN,INC) 215 : * IF (LINCX) INCX.A(INBVA)=VPOCHA(IN,INC) 216 : KS2B.A(IDIAGD)=VPOCHA(IN,INC) 217 : IF (LINCX) INCX.A(IDIAGP)=VPOCHA(IN,INC) 218 : DDLPD.LECT(IDIAGP)=IDIAGD 219 : DDLDP.LECT(IDIAGD)=IDIAGP 220 : * DDLCL.LOGI(INBVA)=.TRUE. 221 : ENDIF 222 : ENDIF 223 : 112 CONTINUE 224 : ENDIF 225 : 11 CONTINUE 226 : SEGDES MSOUPO 227 : SEGDES MELEME 228 : SEGDES MPOVAL 229 : 1 CONTINUE 230 : SEGDES MCHPOI 231 : SEGDES IDMATD 232 : SEGDES IDMATP 233 : SEGSUP MLENTI 234 : SEGDES MINC 235 : IF (IMPR.GT.7) THEN 236 : SEGPRT,KS2B 237 : IF (LINCX) SEGPRT,INCX 238 : ENDIF 239 : IF (LINCX) SEGDES,INCX 240 : C 2) On s'occupe de la matrice : 241 : C pour les no de DDL ou on impose une CL dirichlet, il faut : 242 : C annuler la ligne correspondante sauf la diagonale qui vaut 1 243 : C reporter -colonne*clim sur le second membre 244 : C On parcourt donc toute la matrice 245 : IF (IMPR.GT.6) THEN 246 : WRITE(IOIMP,*) 'Dirichlet sur la matrice' 247 : ENDIF 248 : SEGINI,KMORS=AMORS 249 : SEGINI,KISA=AISA 250 : DO 2 INTTT=1,NTTT 251 : *! IF (DDLCL.LOGI(INTTT)) THEN 252 : IDIAGP=DDLDP.LECT(INTTT) 253 : IF (IDIAGP.NE.0) THEN 254 : C 255 : C CL de dirichlet pour la ligne INTTT : 256 : C On met 1 sur la diagonale et 0 ailleurs, 257 : C le second membre ne change pas 258 : C 259 : FLDIA=.FALSE. 260 : DO 21 ICOL=KMORS.IA(INTTT),(KMORS.IA(INTTT+1)-1) 261 : *! IF (KMORS.JA(ICOL).NE.INTTT) THEN 262 : IF (KMORS.JA(ICOL).NE.IDIAGP) THEN 263 : KISA.A(ICOL)=ZERO 264 : ELSE 265 : KISA.A(ICOL)=ONE 266 : FLDIA=.TRUE. 267 : ENDIF 268 : 21 CONTINUE 269 : IF (.NOT.FLDIA) THEN 270 : *! WRITE(IOIMP,*) 'melim.eso : diag.',INTTT,'inexistante' 271 : WRITE(IOIMP,*) 'melim.eso : diag.',IDIAGP, 272 : $ 'inexistante' 273 : GOTO 9999 274 : ENDIF 275 : ELSE 276 : C 277 : C Ligne INTTT sans condition de dirichlet : 278 : C On annule les termes des colonnes liées à une CL dirichlet, 279 : C et on reporte le produit colonne*CL au second membre 280 : C 281 : TEMP=ZERO 282 : DO 22 ICOL=KMORS.IA(INTTT),(KMORS.IA(INTTT+1)-1) 283 : JNTTT=KMORS.JA(ICOL) 284 : IDIAGD=DDLPD.LECT(JNTTT) 285 : *! IF (DDLCL.LOGI(JNTTT)) THEN 286 : IF (IDIAGD.NE.0) THEN 287 : *! IF (KS2B.A(JNTTT).NE.ZERO) THEN 288 : *! TEMP=TEMP+(KISA.A(ICOL)*KS2B.A(JNTTT)) 289 : TEMP=TEMP+(KISA.A(ICOL)*KS2B.A(IDIAGD)) 290 : *! ENDIF 291 : KISA.A(ICOL)=ZERO 292 : ENDIF 293 : 22 CONTINUE 294 : IF (TEMP.NE.ZERO) THEN 295 : KS2B.A(INTTT)=KS2B.A(INTTT)-TEMP 296 : ENDIF 297 : ENDIF 298 : 2 CONTINUE 299 : SEGDES KISA 300 : SEGDES KMORS 301 : SEGDES KS2B 302 : SEGSUP DDLPD 303 : SEGSUP DDLDP 304 : ENDIF 305 : * 306 : * Terminaison normale 307 : * 308 : IRET=0 309 : RETURN 310 : * 311 : * Format handling 312 : * 313 : 1002 FORMAT(10(1X,1PE11.4)) 314 : * 315 : * Error handling 316 : * 317 : 9999 CONTINUE 318 : IRET=1 319 : WRITE(IOIMP,*) 'An error was detected in melim.eso' 320 : RETURN 321 : * 322 : * End of MELIM 323 : * 324 : END 325 : 326 : 327 : 328 : 329 : 330 : 331 : 332 : 333 : 334 : 335 : 336 : 337 : 338 : 339 : 340 : 341 : 342 : 343 :
© Cast3M 2003 - All rights reserved.
Disclaimer