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