1 : C INELPR SOURCE GOUNAND 06/12/19 21:15:20 5612 2 : SUBROUTINE INELPR(MYLRFS,IMPR,IRET) 3 : IMPLICIT REAL*8 (A-H,O-Z) 4 : IMPLICIT INTEGER (I-N) 5 : C*********************************************************************** 6 : C NOM : INELPR 7 : C PROJET : Noyau linéaire NLIN 8 : C DESCRIPTION : Remplit le segment des éléments de référence 9 : C avec les éléments de référence de dimension 3, 10 : C de forme géométrique prismatique à base triangle. 11 : C 12 : C REFERENCES : 13 : C LANGAGE : ESOPE 14 : C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/LTMF) 15 : C mél : gounand@semt2.smts.cea.fr 16 : C*********************************************************************** 17 : C APPELES : INILRF, INILAG, FILRF, PROLRF, PROBAP, GBAPCO 18 : C APPELE PAR : INLRFS 19 : C*********************************************************************** 20 : C ENTREES : - 21 : C ENTREES/SORTIES : MYLRFS 22 : C SORTIES : - 23 : C CODE RETOUR (IRET) : = 0 si tout s'est bien passé 24 : C*********************************************************************** 25 : C VERSION : v1, 23/03/00, version initiale 26 : C HISTORIQUE : v1, 23/03/00, création 27 : C HISTORIQUE : v2, 10/05/00, modif. du segment ELREF 28 : C HISTORIQUE : 29 : C HISTORIQUE : 30 : C*********************************************************************** 31 : C Prière de PRENDRE LE TEMPS de compléter les commentaires 32 : C en cas de modification de ce sous-programme afin de faciliter 33 : C la maintenance ! 34 : C*********************************************************************** 35 : -INC CCOPTIO 36 : CBEGININCLUDE SELREF 37 : SEGMENT ELREF 38 : CHARACTER*(LNNOM) NOMLRF 39 : CHARACTER*(LNFORM) FORME 40 : CHARACTER*(LNTYPL) TYPEL 41 : CHARACTER*(LNESP) ESPACE 42 : INTEGER DEGRE 43 : REAL*8 XCONOD(NDIMEL,NBNOD) 44 : INTEGER NPQUAF(NBDDL) 45 : INTEGER NUMCMP(NBDDL) 46 : INTEGER QUENOD(NBDDL) 47 : INTEGER ORDDER(NDIMEL,NBDDL) 48 : POINTEUR MBPOLY.POLYNS 49 : ENDSEGMENT 50 : SEGMENT ELREFS 51 : POINTEUR LISEL(0).ELREF 52 : ENDSEGMENT 53 : CENDINCLUDE SELREF 54 : POINTEUR MYLRFS.ELREFS 55 : POINTEUR ELCOUR.ELREF 56 : POINTEUR ELPRO1.ELREF 57 : POINTEUR ELPRO2.ELREF 58 : CBEGININCLUDE SPOLYNO 59 : SEGMENT POLYNO 60 : REAL*8 COEMON(NBMON) 61 : INTEGER EXPMON(NDIML,NBMON) 62 : ENDSEGMENT 63 : SEGMENT POLYNS 64 : POINTEUR LIPOLY(NBPOLY).POLYNO 65 : ENDSEGMENT 66 : CENDINCLUDE SPOLYNO 67 : POINTEUR MYBPOL.POLYNS 68 : POINTEUR MYBPO1.POLYNS 69 : POINTEUR MYBPO2.POLYNS 70 : POINTEUR IZ2.POLYNO 71 : * 72 : INTEGER IMPR,IRET 73 : * Elément de nom : L2D1PR4 74 : REAL*8 UNS6,UNS2,DEUXS3,UNS3,UN,ZERO 75 : PARAMETER (UNS6=1.D0/6.D0) 76 : PARAMETER (UNS2=1.D0/2.D0) 77 : PARAMETER (DEUXS3=2.D0/3.D0) 78 : PARAMETER (UNS3=1.D0/3.D0) 79 : PARAMETER (UN=1.D0) 80 : PARAMETER (ZERO=0.D0) 81 : * 82 : INTEGER INDDL 83 : * 84 : * Executable statements 85 : * 86 : IF (IMPR.GT.1) WRITE(IOIMP,*) 'Entrée dans inelpr' 87 : * 88 : * Elément de nom : L2D0PR1 89 : * Sur un prisme : élément de Lagrange, fonction L2, approximation 90 : * nodale, espace de référence de dimension 3, 1 noeuds, 1 degrés de 91 : * liberté, degré de l'approximation : 0 92 : * 93 : * In INILRF : SEGINI ELCOUR 93 : CALL INILRF('L2D0PR1','PRISME','LAGRANGE','L2', 94 : $ 3,1,1,0, 95 : $ ELCOUR, 96 : $ IMPR,IRET) 97 : IF (IRET.NE.0) GOTO 9999 97 : CALL FILRF('L2D0TR1',MYLRFS,ELPRO1,IMPR,IRET) 98 : IF (IRET.NE.0) GOTO 9999 98 : CALL FILRF('L2D0SE1',MYLRFS,ELPRO2,IMPR,IRET) 99 : IF (IRET.NE.0) GOTO 9999 99 : CALL PROLRF(ELPRO1,ELPRO2,ELCOUR,IMPR,IRET) 100 : IF (IRET.NE.0) GOTO 9999 101 : ELCOUR.NPQUAF(1)=21 102 : ELCOUR.NUMCMP(1)=1 103 : * Initialise la correspondance ddl-noeud+ord.der 103 : CALL INILAG(ELCOUR,IMPR,IRET) 104 : IF (IRET.NE.0) GOTO 9999 105 : SEGACT ELPRO1 106 : MYBPO1=ELPRO1.MBPOLY 107 : SEGDES ELPRO1 108 : SEGACT ELPRO2 109 : MYBPO2=ELPRO2.MBPOLY 110 : SEGDES ELPRO2 111 : * Calcule la base polynômiale produit 111 : CALL PROBAP(MYBPO1,MYBPO2,MYBPOL,IMPR,IRET) 112 : IF (IRET.NE.0) GOTO 9999 113 : ELCOUR.MBPOLY=MYBPOL 114 : SEGDES ELCOUR 115 : MYLRFS.LISEL(**)=ELCOUR 116 : * 117 : * Elément de nom : L2D1PR4 118 : * Sur un prisme : élément de Lagrange, fonction L2, approximation 119 : * nodale, espace de référence de dimension 3, 4 noeuds, 4 degrés de 120 : * liberté, degré de l'approximation : 1 121 : * 122 : * In INILRF : SEGINI ELCOUR 122 : CALL INILRF('L2D1PR4','PRISME','LAGRANGE','L2', 123 : $ 3,4,4,1, 124 : $ ELCOUR, 125 : $ IMPR,IRET) 126 : IF (IRET.NE.0) GOTO 9999 127 : ELCOUR.XCONOD(1,1)= UNS6 128 : ELCOUR.XCONOD(2,1)= UNS6 129 : ELCOUR.XCONOD(3,1)=-UNS2 130 : ELCOUR.XCONOD(1,2)= DEUXS3 131 : ELCOUR.XCONOD(2,2)= UNS6 132 : ELCOUR.XCONOD(3,2)= UNS2 133 : ELCOUR.XCONOD(1,3)= UNS6 134 : ELCOUR.XCONOD(2,3)= DEUXS3 135 : ELCOUR.XCONOD(3,3)=-UNS2 136 : ELCOUR.XCONOD(1,4)= UNS3 137 : ELCOUR.XCONOD(2,4)= UNS3 138 : ELCOUR.XCONOD(3,4)= UNS2 139 : * Les d.d.l. sont au noeud 21 140 : DO 199 INDDL=1,4 141 : ELCOUR.NPQUAF(INDDL)=21 142 : ELCOUR.NUMCMP(INDDL)=INDDL 143 : 199 CONTINUE 144 : * Initialise la correspondance ddl-noeud+ord.der 144 : CALL INILAG(ELCOUR,IMPR,IRET) 145 : IF (IRET.NE.0) GOTO 9999 146 : * Génère une base polynômiale complète (dimension 3, degré 1) 146 : CALL GBAPCO(3,1,MYBPOL,IMPR,IRET) 147 : IF (IRET.NE.0) GOTO 9999 148 : ELCOUR.MBPOLY=MYBPOL 149 : SEGDES ELCOUR 150 : MYLRFS.LISEL(**)=ELCOUR 151 : * 152 : * Elément de nom : H1D1PR6 153 : * Sur un prisme : élément de Lagrange, fonction H1, approximation 154 : * nodale, espace de référence de dimension 3, 6 noeuds, 6 degrés de 155 : * liberté, degré de l'approximation : 1 156 : * 157 : * In INILRF : SEGINI ELCOUR 157 : CALL INILRF('H1D1PR6','PRISME','LAGRANGE','H1', 158 : $ 3,6,6,1, 159 : $ ELCOUR, 160 : $ IMPR,IRET) 161 : IF (IRET.NE.0) GOTO 9999 161 : CALL FILRF('H1D1TR3',MYLRFS,ELPRO1,IMPR,IRET) 162 : IF (IRET.NE.0) GOTO 9999 162 : CALL FILRF('H1D1SE2',MYLRFS,ELPRO2,IMPR,IRET) 163 : IF (IRET.NE.0) GOTO 9999 163 : CALL PROLRF(ELPRO1,ELPRO2,ELCOUR,IMPR,IRET) 164 : IF (IRET.NE.0) GOTO 9999 165 : * Les d.d.l. sont aux noeuds 1,3,5... 166 : DO 201 INDDL=1,3 167 : ELCOUR.NPQUAF(INDDL)=(2*INDDL)-1 168 : ELCOUR.NUMCMP(INDDL)=1 169 : 201 CONTINUE 170 : * ... et 10,12,14. 171 : DO 203 INDDL=4,6 172 : ELCOUR.NPQUAF(INDDL)=(2*(INDDL-4))+10 173 : ELCOUR.NUMCMP(INDDL)=1 174 : 203 CONTINUE 175 : * Initialise la correspondance ddl-noeud+ord.der 175 : CALL INILAG(ELCOUR,IMPR,IRET) 176 : IF (IRET.NE.0) GOTO 9999 177 : SEGACT ELPRO1 178 : MYBPO1=ELPRO1.MBPOLY 179 : SEGDES ELPRO1 180 : SEGACT ELPRO2 181 : MYBPO2=ELPRO2.MBPOLY 182 : SEGDES ELPRO2 183 : * Calcule la base polynômiale produit 183 : CALL PROBAP(MYBPO1,MYBPO2,MYBPOL,IMPR,IRET) 184 : IF (IRET.NE.0) GOTO 9999 185 : ELCOUR.MBPOLY=MYBPOL 186 : SEGDES ELCOUR 187 : MYLRFS.LISEL(**)=ELCOUR 188 : * 189 : * Elément de nom : CRD1PR5 190 : * Sur un prisme : élément de Lagrange, fonction continue au centre des 191 : * faces, approximation nodale, espace de référence de dimension 3, 5 192 : * noeuds, 5 degrés de 193 : * liberté, degré de l'approximation : 1 194 : * 195 : * In INILRF : SEGINI ELCOUR 195 : CALL INILRF('CRD1PR5','PRISME','LAGRANGE','HFAC', 196 : $ 3,5,5,1, 197 : $ ELCOUR, 198 : $ IMPR,IRET) 199 : IF (IRET.NE.0) GOTO 9999 200 : ELCOUR.XCONOD(1,1)=UNS2 201 : ELCOUR.XCONOD(2,1)=ZERO 202 : ELCOUR.XCONOD(3,1)=ZERO 203 : ELCOUR.XCONOD(1,2)=UNS2 204 : ELCOUR.XCONOD(2,2)=UNS2 205 : ELCOUR.XCONOD(3,2)=ZERO 206 : ELCOUR.XCONOD(1,3)=ZERO 207 : ELCOUR.XCONOD(2,3)=UNS2 208 : ELCOUR.XCONOD(3,3)=ZERO 209 : ELCOUR.XCONOD(1,4)=UNS3 210 : ELCOUR.XCONOD(2,4)=UNS3 211 : ELCOUR.XCONOD(3,4)=-UN 212 : ELCOUR.XCONOD(1,5)=UNS3 213 : ELCOUR.XCONOD(2,5)=UNS3 214 : ELCOUR.XCONOD(3,5)=+UN 215 : * Les d.d.l. sont aux noeuds 16,17,18,19,20 216 : DO INDDL=1,5 217 : ELCOUR.NPQUAF(INDDL)=INDDL+15 218 : ELCOUR.NUMCMP(INDDL)=1 219 : ENDDO 220 : * Initialise la correspondance ddl-noeud+ord.der 220 : CALL INILAG(ELCOUR,IMPR,IRET) 221 : IF (IRET.NE.0) GOTO 9999 222 : * Génère une base polynômiale complète (dimension 3, degré 1) 222 : CALL GBAPCO(3,1,MYBPOL,IMPR,IRET) 223 : IF (IRET.NE.0) GOTO 9999 224 : * Il faut rajouter le monôme z^2 225 : SEGACT MYBPOL*MOD 226 : NDIML=3 227 : NBMON=1 228 : SEGINI IZ2 229 : IZ2.COEMON(1)=UN 230 : IZ2.EXPMON(3,1)=2 231 : SEGDES IZ2 232 : MYBPOL.LIPOLY(**)=IZ2 233 : SEGDES MYBPOL 234 : ELCOUR.MBPOLY=MYBPOL 235 : SEGDES ELCOUR 236 : MYLRFS.LISEL(**)=ELCOUR 237 : * 238 : * Elément de nom : H1D2PR15 239 : * Sur un prisme : élément de Lagrange incomplet (Serendip), 240 : * fonction H1, approximation 241 : * nodale, espace de référence de dimension 3, 15 noeuds, 15 degrés de 242 : * liberté, degré de l'approximation : 2 243 : * 244 : * In INILRF : SEGINI ELCOUR 244 : CALL INILRF('H1D2PR15','PRISME','LAGRANGE','H1', 245 : $ 3,15,15,2, 246 : $ ELCOUR, 247 : $ IMPR,IRET) 248 : IF (IRET.NE.0) GOTO 9999 248 : CALL FILRF('H1D2TR6',MYLRFS,ELPRO1,IMPR,IRET) 249 : IF (IRET.NE.0) GOTO 9999 249 : CALL FILRF('H1D1TR3',MYLRFS,ELPRO2,IMPR,IRET) 250 : IF (IRET.NE.0) GOTO 9999 251 : C Inutile 252 : C ELCOUR.XCONOD(1,1)=ZERO 253 : C ELCOUR.XCONOD(2,1)=ZERO 254 : C ELCOUR.XCONOD(3,1)=ZERO 255 : * Les d.d.l. sont aux noeuds 1,...,15 256 : DO 217 INDDL=1,15 257 : ELCOUR.NPQUAF(INDDL)=INDDL 258 : ELCOUR.NUMCMP(INDDL)=1 259 : 217 CONTINUE 260 : * Initialise la correspondance ddl-noeud+ord.der 260 : CALL INILAG(ELCOUR,IMPR,IRET) 261 : IF (IRET.NE.0) GOTO 9999 262 : * Pas de base polynômiale (on recopie l'élément de castem) 263 : ELCOUR.MBPOLY=0 264 : SEGDES ELCOUR 265 : MYLRFS.LISEL(**)=ELCOUR 266 : * 267 : * Elément de nom : H1D2PR18 268 : * Sur un prisme : élément de Lagrange, fonction H1, approximation 269 : * nodale, espace de référence de dimension 3, 18 noeuds, 18 degrés de 270 : * liberté, degré de l'approximation : 2 271 : * 272 : * In INILRF : SEGINI ELCOUR 272 : CALL INILRF('H1D2PR18','PRISME','LAGRANGE','H1', 273 : $ 3,18,18,2, 274 : $ ELCOUR, 275 : $ IMPR,IRET) 276 : IF (IRET.NE.0) GOTO 9999 276 : CALL FILRF('H1D2TR6',MYLRFS,ELPRO1,IMPR,IRET) 277 : IF (IRET.NE.0) GOTO 9999 277 : CALL FILRF('H1D2SE3',MYLRFS,ELPRO2,IMPR,IRET) 278 : IF (IRET.NE.0) GOTO 9999 278 : CALL PROLRF(ELPRO1,ELPRO2,ELCOUR,IMPR,IRET) 279 : IF (IRET.NE.0) GOTO 9999 280 : * Les ieme d.d.l sont aux noeuds j 281 : DO IDDL=1,6 282 : ELCOUR.NPQUAF(IDDL)=IDDL 283 : ENDDO 284 : ELCOUR.NPQUAF( 7)= 7 285 : ELCOUR.NPQUAF( 8)=16 286 : ELCOUR.NPQUAF( 9)= 8 287 : ELCOUR.NPQUAF(10)=17 288 : ELCOUR.NPQUAF(11)= 9 289 : ELCOUR.NPQUAF(12)=18 290 : DO IDDL=13,18 291 : ELCOUR.NPQUAF(IDDL)=IDDL-3 292 : ENDDO 293 : DO IDDL=1,18 294 : ELCOUR.NUMCMP(IDDL)=1 295 : ENDDO 296 : * Initialise la correspondance ddl-noeud+ord.der 296 : CALL INILAG(ELCOUR,IMPR,IRET) 297 : IF (IRET.NE.0) GOTO 9999 298 : SEGACT ELPRO1 299 : MYBPO1=ELPRO1.MBPOLY 300 : SEGDES ELPRO1 301 : SEGACT ELPRO2 302 : MYBPO2=ELPRO2.MBPOLY 303 : SEGDES ELPRO2 304 : * Calcule la base polynômiale produit 304 : CALL PROBAP(MYBPO1,MYBPO2,MYBPOL,IMPR,IRET) 305 : IF (IRET.NE.0) GOTO 9999 306 : ELCOUR.MBPOLY=MYBPOL 307 : SEGDES ELCOUR 308 : MYLRFS.LISEL(**)=ELCOUR 309 : * 310 : * Elément de nom : H1D2PR21 311 : * Sur un prisme : élément de Lagrange, fonction H1, approximation 312 : * nodale, espace de référence de dimension 3, 21 noeuds, 21 degrés de 313 : * liberté, degré de l'approximation : 2 314 : * 315 : * In INILRF : SEGINI ELCOUR 315 : CALL INILRF('H1D2PR21','PRISME','LAGRANGE','H1', 316 : $ 3,21,21,2, 317 : $ ELCOUR, 318 : $ IMPR,IRET) 319 : IF (IRET.NE.0) GOTO 9999 319 : CALL FILRF('H1D2TR7',MYLRFS,ELPRO1,IMPR,IRET) 320 : IF (IRET.NE.0) GOTO 9999 320 : CALL FILRF('H1D2SE3',MYLRFS,ELPRO2,IMPR,IRET) 321 : IF (IRET.NE.0) GOTO 9999 321 : CALL PROLRF(ELPRO1,ELPRO2,ELCOUR,IMPR,IRET) 322 : IF (IRET.NE.0) GOTO 9999 323 : * Les ieme d.d.l sont aux noeuds j 324 : DO IDDL=1,6 325 : ELCOUR.NPQUAF(IDDL)=IDDL 326 : ENDDO 327 : ELCOUR.NPQUAF( 7)=19 328 : ELCOUR.NPQUAF( 8)= 7 329 : ELCOUR.NPQUAF( 9)=16 330 : ELCOUR.NPQUAF(10)= 8 331 : ELCOUR.NPQUAF(11)=17 332 : ELCOUR.NPQUAF(12)= 9 333 : ELCOUR.NPQUAF(13)=18 334 : ELCOUR.NPQUAF(14)=21 335 : DO IDDL=15,20 336 : ELCOUR.NPQUAF(IDDL)=IDDL-5 337 : ENDDO 338 : ELCOUR.NPQUAF(21)=20 339 : DO IDDL=1,21 340 : ELCOUR.NUMCMP(IDDL)=1 341 : ENDDO 342 : * Initialise la correspondance ddl-noeud+ord.der 342 : CALL INILAG(ELCOUR,IMPR,IRET) 343 : IF (IRET.NE.0) GOTO 9999 344 : SEGACT ELPRO1 345 : MYBPO1=ELPRO1.MBPOLY 346 : SEGDES ELPRO1 347 : SEGACT ELPRO2 348 : MYBPO2=ELPRO2.MBPOLY 349 : SEGDES ELPRO2 350 : * Calcule la base polynômiale produit 350 : CALL PROBAP(MYBPO1,MYBPO2,MYBPOL,IMPR,IRET) 351 : IF (IRET.NE.0) GOTO 9999 352 : ELCOUR.MBPOLY=MYBPOL 353 : SEGDES ELCOUR 354 : MYLRFS.LISEL(**)=ELCOUR 355 : * 356 : * Normal termination 357 : * 358 : IRET=0 359 : RETURN 360 : * 361 : * Format handling 362 : * 363 : * 364 : * Error handling 365 : * 366 : 9999 CONTINUE 367 : IRET=1 368 : WRITE(IOIMP,*) 'An error was detected in subroutine inelpr' 369 : RETURN 370 : * 371 : * End of subroutine INELPR 372 : * 373 : END 374 : 375 : 376 : 377 :
© Cast3M 2003 - All rights reserved.
Disclaimer