1 : C TGRAD SOURCE CHAT 09/10/09 21:25:20 6519 2 : SUBROUTINE TGRAD (IPTAB1) 3 : C 4 : C*********************************************************************** 5 : C 6 : C 7 : C FONCTION: 8 : C --------- 9 : C 10 : C en "0D/1D" sur des éléments de type POINT. 11 : C 12 : C 13 : C ENTREE : 14 : C -------- 15 : C 16 : C IPTAB1 : TABLE de soustype 'OPER_0D' contenant les indices 17 : C suivants : 18 : C 19 : C TAB1 . 'GEOINF' : TABLE des informations géométriques de soustype 20 : C 'GEOINF' (type ENTIER). 21 : C TAB1 . 'INCO' : TABLE de soustype 'INCO' contenant l'ensemble 22 : C des champs à l'itération précédant l'itération 23 : C courante (type ENTIER). 24 : C TAB1 . 'POTENTIA' : Valeurs du champ potentiel (gravité * cote) 25 : C (type CHPOINT) (indice facultatif). 26 : C Cet indice est utilisé pour calculer p.e le 27 : C gradient de rho * gravité * cote. 28 : C (doit avoir le même support que l'inconnue 29 : C duale). 30 : C TAB1 . 'DUAL' : Nom de l'inconnue duale (doit être un indice de 31 : C la table de soustype 'INCO') (type MOT). 32 : C TAB1 . 'PRIMAL' : Nom de l'inconnue primale (doit être un indice 33 : C de la table de soustype 'INCO') (type MOT). 34 : C 35 : C 36 : C RESULTAT : 37 : C ---------- 38 : C 39 : C TAB1 . 'LHS' : Matrice élémentaire associée à l'opération 40 : C (type RIGIDITE). 41 : C 42 : C 43 : C AUTEUR, DATE DE CREATION: 44 : C ------------------------- 45 : C 46 : C Laurent DADA décembre 1996 47 : C 48 : C 49 : C LANGAGE: 50 : C -------- 51 : C 52 : C ESOPE + FORTRAN77 53 : C 54 : C*********************************************************************** 55 : C 56 : IMPLICIT INTEGER(I-N) 57 : IMPLICIT REAL*8 (A-H,O-Z) 58 : C 59 : -INC CCOPTIO 60 : -INC CCGEOME 61 : -INC SMCOORD 62 : -INC SMTABLE 63 : POINTEUR MTABG.MTABLE,MTABS.MTABLE,IPTABI.MTABLE 64 : -INC SMCHPOI 65 : -INC SMELEME 66 : POINTEUR IPTD1.MELEME,IPTP1.MELEME,IPJUCE.MELEME 67 : -INC SMRIGID 68 : C 69 : SEGMENT REDI 70 : INTEGER ORDR1(NNGOT) 71 : INTEGER ORDR2(NNGOT) 72 : ENDSEGMENT 73 : C 74 : CHARACTER*8 TYPE,MOTI,MOT1,NOMPR1,NOMDU1,NOSUP1,NOSUD1 75 : CHARACTER*7 NAMT1 76 : LOGICAL LPOT 77 : C 78 : C Lecture de la table GEOINF de la table OPER_0D 79 : C 80 : TYPE = 'TABLE ' 80 : CALL ACMO (IPTAB1,'GEOINF',TYPE,IPTABG) 81 : IF (IERR.NE.0) RETURN 82 : C 83 : MOTI = 'SOUSTYPE' 83 : CALL ACMM (IPTABG,MOTI,MOT1) 84 : IF (IERR.NE.0) RETURN 85 : IF (MOT1(1:6).NE.'GEOINF') THEN 86 : MOTERR(1:8) = 'GEOINF ' 87 : MOTERR(9:16) = 'GEOINF ' 87 : CALL ERREUR (790) 88 : RETURN 89 : ENDIF 90 : C 91 : C Lecture de la table INCO dans la table OPER_0D 92 : C 93 : TYPE = 'TABLE ' 93 : CALL ACMO (IPTAB1,'INCO',TYPE,IPTABI) 94 : IF (IERR.NE.0) RETURN 95 : C 96 : MOTI = 'SOUSTYPE' 96 : CALL ACMM (IPTABI,MOTI,MOT1) 97 : IF (IERR.NE.0) RETURN 98 : IF (MOT1(1:4).NE.'INCO') THEN 99 : MOTERR(1:8) = 'INCO ' 100 : MOTERR(9:16) = 'INCO ' 100 : CALL ERREUR (790) 101 : RETURN 102 : ENDIF 103 : C 104 : C Lecture de la table SUPPORT dans la table INCO 105 : C 106 : TYPE = 'TABLE ' 106 : CALL ACMO (IPTABI,'SUPPORT',TYPE,IPTABS) 107 : IF (IERR.NE.0) RETURN 108 : C 109 : C Lecture du MAILLAGE des connectivités 'JUNCEL' de la table GEOINF 110 : C Arrêt si les éléments ne sont pas des SEG3 111 : C 112 : TYPE = 'MAILLAGE' 112 : CALL ACMO (IPTABG,'JUNCEL',TYPE,IPJUCE) 113 : IF (IERR.NE.0) RETURN 114 : SEGACT IPJUCE 115 : IF ((IPJUCE.ITYPEL).NE.3) THEN 116 : MOTERR(1:8) = 'JUNCEL ' 117 : MOTERR(9:16) = 'MAILLAGE' 117 : CALL ERREUR (787) 118 : SEGDES IPJUCE 119 : RETURN 120 : ENDIF 121 : NBEJC1 = IPJUCE.NUM(/2) 122 : C 123 : C Lecture du nom de l'inconnue PRIMAL 124 : C 125 : TYPE = ' ' 125 : CALL ACMO (IPTAB1,'PRIMAL',TYPE,IPR1) 126 : IF (IERR.NE.0) RETURN 127 : IF (TYPE.EQ.'MOT ') THEN 127 : CALL ACMM (IPTAB1,'PRIMAL',NOMPR1) 128 : IF (IERR.NE.0) RETURN 129 : ELSE 130 : MOTERR(1:8) = 'PRIMAL ' 131 : MOTERR(9:16) = TYPE 131 : CALL ERREUR (787) 132 : RETURN 133 : ENDIF 134 : C 135 : C Lecture du nom de l'inconnue DUAL 136 : C 137 : TYPE = ' ' 137 : CALL ACMO (IPTAB1,'DUAL',TYPE,IDU1) 138 : IF (IERR.NE.0) RETURN 139 : IF (TYPE.EQ.'MOT ') THEN 139 : CALL ACMM (IPTAB1,'DUAL',NOMDU1) 140 : IF (IERR.NE.0) RETURN 141 : ELSE 142 : MOTERR(1:8) = 'DUAL ' 143 : MOTERR(9:16) = TYPE 143 : CALL ERREUR (787) 144 : RETURN 145 : ENDIF 146 : C 147 : C Lecture éventuelle du champ POTENTIA 148 : C 149 : LPOT = .FALSE. 150 : TYPE = ' ' 150 : CALL ACMO (IPTAB1,'POTENTIA',TYPE,IPCHP1) 151 : IF (IERR.NE.0) RETURN 152 : IF (TYPE.EQ.'CHPOINT ') THEN 153 : LPOT = .TRUE. 154 : C Contrôle du CHPOINT 155 : MCHPOI = IPCHP1 156 : SEGACT MCHPOI 157 : MSOUPO = IPCHP(1) 158 : SEGDES MCHPOI 159 : SEGACT MSOUPO 160 : NC = NOHARM(/1) 161 : IF (NC.NE.1) THEN 162 : MOTERR(1:8) = 'POTENTIA' 163 : MOTERR(9:16) = 'CHPOINT ' 163 : CALL ERREUR (784) 164 : SEGDES MSOUPO 165 : RETURN 166 : ENDIF 167 : IPT1 = IGEOC 168 : MPOVA1 = IPOVAL 169 : SEGDES MSOUPO 170 : ENDIF 171 : C 172 : C Lecture du nom du support de l'inconnue PRIMAL 173 : C Lecture du MAILLAGE de l'inconnue PRIMAL 174 : C 175 : TYPE = ' ' 175 : CALL ACMO (IPTABS,NOMPR1,TYPE,ISUP1) 176 : IF (IERR.NE.0) RETURN 177 : IF (TYPE.EQ.'MOT ') THEN 177 : CALL ACMM (IPTABS,NOMPR1,NOSUP1) 178 : IF (IERR.NE.0) RETURN 179 : ENDIF 180 : TYPE = 'MAILLAGE' 180 : CALL ACMO (IPTABG,NOSUP1,TYPE,IPTP1) 181 : IF (IERR.NE.0) RETURN 182 : C 183 : C Lecture du nom du support de l'inconnue DUAL 184 : C Arrêt si différent de 'JUNCTION' 185 : C Lecture du MAILLAGE de l'inconnue DUAL 186 : C 187 : TYPE = ' ' 187 : CALL ACMO (IPTABS,NOMDU1,TYPE,ISUD1) 188 : IF (TYPE.EQ.'MOT ') THEN 188 : CALL ACMM (IPTABS,NOMDU1,NOSUD1) 189 : IF (IERR.NE.0) RETURN 190 : IF (NOSUD1.NE.'JUNCTION') THEN 191 : MOTERR(1:8) = 'DUAL ' 192 : MOTERR(9:16) = 'CHPOINT ' 192 : CALL ERREUR (788) 193 : RETURN 194 : ENDIF 195 : ENDIF 196 : TYPE = 'MAILLAGE' 196 : CALL ACMO (IPTABG,NOSUD1,TYPE,IPTD1) 197 : IF (IERR.NE.0) RETURN 198 : C 199 : C Création du support géométrique pour la RIGIDITE 200 : C (maillage de type SUPER-ELEMENT). 201 : C 202 : NNGOT = XCOOR(/1)/(IDIM+1) 203 : SEGINI REDI 204 : C 205 : C On fusionne les maillages de POI1 des supports des inconnues PRIMAL 206 : C et DUAL en un maillage de type SUPER-ELEMENT. 207 : C 208 : SEGACT IPTD1 208 : IF (IPTD1.ITYPEL.NE.1) CALL CHANGE (IPTD1,1) 209 : SEGACT IPTD1 210 : SEGACT IPTP1 210 : IF (IPTP1.ITYPEL.NE.1) CALL CHANGE (IPTP1,1) 211 : SEGACT IPTP1 212 : NBNNP1 = IPTP1.NUM(/2) 213 : NBNND1 = IPTD1.NUM(/2) 214 : NBNN = NBNNP1 + NBNND1 215 : NBSOUS = 0 216 : NBREF = 0 217 : NBELEM = 1 218 : SEGINI MELEME 219 : ICOLOR(1) = IDCOUL 220 : ITYPEL = 28 221 : DO 50 I50=1,NBNNP1 222 : NUM(I50,1) = IPTP1.NUM(1,I50) 223 : ORDR2(NUM(I50,1)) = I50 224 : 50 CONTINUE 225 : DO 60 I60=1,NBNND1 226 : NUM(I60+NBNNP1,1) = IPTD1.NUM(1,I60) 227 : ORDR1(NUM(I60+NBNNP1,1)) = I60 228 : 60 CONTINUE 229 : C 230 : C Création de la RIGIDITE 231 : C 232 : NRIGE = 7 233 : NRIGEL = 1 234 : SEGINI MRIGID 235 : C 236 : MTYMAT = 'RIGIDITE' 237 : IFORIG = IFOMOD 238 : ICHOLE = 0 239 : IMGEO1 = 0 240 : IMGEO2 = 0 241 : ISUPEQ = 0 242 : COERIG(1) = 1.D0 243 : IRIGEL(1,1) = MELEME 244 : IRIGEL(2,1) = 0 245 : IRIGEL(5,1) = NIFOUR 246 : IRIGEL(6,1) = 0 247 : IRIGEL(7,1) = 2 248 : C 249 : SEGDES MELEME 250 : C 251 : C Remplissage du descripteur de l'objet RIGIDITE 252 : C 253 : NLIGRP = NBNNP1 254 : NLIGRD = NBNND1 255 : SEGINI DESCR 256 : IRIGEL(3,1) = DESCR 257 : DO 10 I10=1,NBNNP1 258 : NOELEP(I10) = I10 259 : LISINC(I10) = NOMPR1 260 : 10 CONTINUE 261 : DO 11 I11=1,NBNND1 262 : NOELED(I11) = I11+NBNNP1 263 : LISDUA(I11) = NOMDU1 264 : 11 CONTINUE 265 : C 266 : SEGDES DESCR 267 : C 268 : NELRIG = 1 269 : * SEGINI IMATRI 270 : SEGINI XMATRI 271 : C 272 : * IMATTT(1) = XMATRI 273 : IRIGEL(4,1) = xMATRI 274 : * SEGDES IMATRI 275 : C 276 : C Remplissage de la matrice élémentaire 277 : C 278 : C mise à 0.d0 de la matrice élémentaire 279 : DO 70 I70=1,NBNND1 280 : DO 701 I701=1,NBNNP1 281 : RE(I70,I701,1) = 0.D0 282 : 701 CONTINUE 283 : 70 CONTINUE 284 : C 285 : C avec prise en compte du terme POTENTIA 286 : C 287 : IF (LPOT) THEN 288 : C activation du maillage et des valeurs du CHPOINT POTENTIA 289 : SEGACT IPT1 290 : SEGACT MPOVA1 291 : C petit controle du support du champ POTENTIA 292 : NBEL1 = IPT1.NUM(/2) 293 : IF (NBEL1.NE.(NBNNP1+NBNND1)) THEN 294 : MOTERR(1:8) = 'POTENTIA' 295 : MOTERR(9:16) = 'CHPOINT ' 295 : CALL ERREUR (788) 296 : SEGDES IPT1 297 : SEGDES MPOVA1 298 : SEGDES IPJUCE 299 : SEGSUP MELEME 300 : SEGSUP XMATRI 301 : SEGSUP MRIGID 302 : SEGDES IPTP1 303 : SEGDES IPTD1 304 : SEGSUP REDI 305 : RETURN 306 : ENDIF 307 : DO 72 I72=1,NBEJC1 308 : NPT1 = IPJUCE.NUM(1,I72) 309 : NPTF1 = IPJUCE.NUM(2,I72) 310 : NPT2 = IPJUCE.NUM(3,I72) 311 : C récupération des valeurs du champ POTENTIA aux points NPT1,NPT2,NPTF1 312 : C balayage sur les points du CHPOINT POTENTIA 313 : XGZP1 = 0.D0 314 : XGZD1 = 0.D0 315 : XGZP2 = 0.D0 316 : DO 723 I723=1,NBEL1 317 : IF (NPT1.EQ.(IPT1.NUM(1,I723))) XGZP1 = MPOVA1.VPOCHA(I723,1) 318 : IF (NPT2.EQ.(IPT1.NUM(1,I723))) XGZP2 = MPOVA1.VPOCHA(I723,1) 319 : IF (NPTF1.EQ.(IPT1.NUM(1,I723))) XGZD1 = MPOVA1.VPOCHA(I723,1) 320 : 723 CONTINUE 321 : C positions dans la matrice élémentaire 322 : IPRI1 = ORDR2(NPT1) 323 : IDUA1 = ORDR1(NPTF1) 324 : IPRI2 = ORDR2(NPT2) 325 : C remplissage de la matrice 326 : RE(IDUA1,IPRI1,1) = -1.D0 * (XGZP1 - XGZD1) 327 : RE(IDUA1,IPRI2,1) = (XGZP2 - XGZD1) 328 : 72 CONTINUE 329 : SEGDES IPTP1 330 : SEGDES IPTD1 331 : C 332 : C 333 : ELSE 334 : C 335 : C 336 : DO 82 I82=1,NBEJC1 337 : NPT1 = IPJUCE.NUM(1,I82) 338 : NPTF1 = IPJUCE.NUM(2,I82) 339 : NPT2 = IPJUCE.NUM(3,I82) 340 : C positions dans la matrice élémentaire 341 : IPRI1 = ORDR2(NPT1) 342 : IDUA1 = ORDR1(NPTF1) 343 : IPRI2 = ORDR2(NPT2) 344 : C remplissage de la matrice 345 : RE(IDUA1,IPRI1,1) = -1.D0 346 : RE(IDUA1,IPRI2,1) = 1.D0 347 : 82 CONTINUE 348 : ENDIF 349 : C 350 : SEGDES IPJUCE 351 : SEGDES XMATRI 352 : SEGSUP REDI 353 : C 354 : C Ecriture du résultat 355 : C 356 : TYPE = 'RIGIDITE' 356 : CALL ECMO (IPTAB1,'LHS',TYPE,MRIGID) 357 : IF (IERR.NE.0) RETURN 358 : C 359 : SEGDES MRIGID 360 : C 361 : END 362 : 363 : 364 :
© Cast3M 2003 - All rights reserved.
Disclaimer