1 : C KDOM10 SOURCE CHAT 05/01/13 00:53:29 5004 2 : SUBROUTINE KDOM10(MTAB) 3 : C 4 : C************************************************************************ 5 : C 6 : C PROJET : CASTEM 2000 7 : C 8 : C NOM : KDOM10 9 : C 9 : C DESCRIPTION : Subroutine called by KDOM1 in the case of EULER 10 : C model 11 : C We fill the domain table MTAB 12 : C 13 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec estensions CISI) 14 : C 15 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF 16 : C 17 : C************************************************************************ 18 : C 19 : C 20 : C**** Variables de COOPTIO 21 : C 22 : IMPLICIT INTEGER(I-N) 23 : INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 24 : & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 25 : & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 26 : & ,IECHO, IIMPI, IOSPI 27 : & ,IDIM, IFICLE, IPREFI 28 : C & ,MCOORD 29 : & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 30 : & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 31 : & ,NORINC,NORVAL,NORIND,NORVAD 32 : & ,NUCROU, IPSAUV, IREFOR, ISAFOR 33 : C 34 : INTEGER MTAB, JG, IELEM, ISOUS, NBELEM, NBNN, NBREF, NBSOUS, NBS 35 : & ,NBE0, IGEOM, JGM, JGN, NN1, NN2, NN3, NN4, NN5, NN6, NTP 36 : & ,LAST, LISFAC(6), INOE, NNOEU, NFAC, NNOEUT, LISFAT(4), NFACT 37 : & ,LASTT,IFAC, NSOM, LASTS, NNS, LISSOM(8),LIFAC(5,6) 38 : & ,LIFACT(4,4), NBSFP, MCHPSU, MCHPNO, MCHPMR, MCHDIA 39 : REAl*8 VOL, CEL(3), VOL1 40 : LOGICAL LOGCOM 41 : CHARACTER*8 TYPI 42 : -INC CCOPTIO 43 : -INC SMELEME 44 : -INC SMLENTI 45 : -INC SMCHPOI 46 : -INC SMLMOTS 47 : -INC SMCOORD 48 : POINTEUR MELCEN.MELEME, MELMAI.MELEME, MLREST.MLENTI 49 : & , MLRES.MLENTI, MELTFA.MELEME, MLRESS.MLENTI 50 : & , MELSOM.MELEME, MELFAC.MELEME, MELFAL.MELEME, MLEFAC.MLENTI 51 : & , MLETOF.MLENTI, IPTT.MELEME, IPTQ.MELEME, MELFAP.MELEME 52 : C 53 : C Elements allowed are 'SEG3' 'TRI7' 'QUA9' 'TE15' 'PY19' 'PR21' 'CU27' 54 : C ITYPEL 3 7 11 35 36 34 33 55 : C 56 : TYPI='MAILLAGE' 56 : CALL ACMO(MTAB,'QUAF',TYPI,MELEME) 57 : IF(IERR.NE.0)GOTO 9999 58 : C 59 : C 2D 3D 60 : C**** 'MAILLAGE' complex 1D/2D elts 2D/3D elts 61 : C 'CENTRE' 'POI1' (ITYPEL = 1) 1D/2D elts 2D/3D elts 62 : C 'SOMMET' 'POI1' (ITYPEL = 1) 2D elts 3D elts 63 : C 'FACE' 'POI1' (ITYPEL = 1) 2D elts 3D elts 64 : C 'FACEL' 'SEG3' (ITYPEL = 3) 2D elts 3D elts 65 : C 'FACEP' 'SEG3' (ITYPEL = 3) in 2D 2D elts 3D elts 66 : C complex in 3D 67 : C 68 : C 'ELTFA' complex 2D elts 3D elts 69 : C 70 : SEGACT MELEME 71 : NBS=MELEME.LISOUS(/1) 72 : IF(NBS .EQ. 0) NBS=1 73 : JG=NBS 74 : C MLENTI contains ITYPEL 75 : C MLENT1 contains NBELEM 76 : SEGINI MLENTI 77 : SEGINI MLENT1 78 : DO ISOUS=1,NBS,1 79 : IF(NBS .NE. 1)THEN 80 : IPT1=MELEME.LISOUS(ISOUS) 81 : SEGACT IPT1 82 : ELSE 83 : IPT1=MELEME 84 : ENDIF 85 : C We have already checked the elements available in KDOM2 86 : C We have already checked that they are referred just ONCE 87 : C (each MELEME.LISOUS(ISOUS) has different ITYPEL) 88 : C 89 : MLENTI.LECT(ISOUS)=IPT1.ITYPEL 90 : MLENT1.LECT(ISOUS)=IPT1.NUM(/2) 91 : ENDDO 92 : C 93 : C**** Compatibility check 94 : C 95 : C We do not allow 1D/2D meshes and 2D/3D meshes 96 : C 97 : IF(IDIM.EQ.2)THEN 98 : IF(MLENTI.LECT(1) .EQ. 3)THEN 99 : LOGCOM=.FALSE. 100 : IF(NBS .NE. 1)THEN 101 : WRITE(IOIMP,*) 'Euler model' 101 : CALL ERREUR(21) 102 : GOTO 9999 103 : ENDIF 104 : ELSE 105 : LOGCOM=.TRUE. 106 : DO ISOUS=2,NBS,1 107 : IF(MLENTI.LECT(ISOUS) .EQ. 3)THEN 108 : WRITE(IOIMP,*) 'Euler model' 109 : WRITE(IOIMP,*) 'No mixing between 1D and 2D' 109 : CALL ERREUR(21) 110 : GOTO 9999 111 : ENDIF 112 : ENDDO 113 : ENDIF 114 : ELSE 115 : C 116 : C****** 3D. 117 : C No 1D possibility 118 : C 119 : DO ISOUS=1,NBS,1 120 : IF(MLENTI.LECT(ISOUS) .EQ. 3)THEN 121 : WRITE(IOIMP,*) 'Euler model' 122 : WRITE(IOIMP,*) 'No 1D meshes in 3D' 122 : CALL ERREUR(21) 123 : GOTO 9999 124 : ENDIF 125 : ENDDO 126 : C 127 : IF((MLENTI.LECT(1) .EQ. 7) .OR. (MLENTI.LECT(1) .EQ. 11))THEN 128 : LOGCOM=.FALSE. 129 : DO ISOUS=1,NBS,1 130 : IF((MLENTI.LECT(ISOUS) .NE. 7) .AND. 131 : & (MLENTI.LECT(ISOUS) .NE. 11))THEN 132 : WRITE(IOIMP,*) 'Euler model' 133 : WRITE(IOIMP,*) 'No mixing between 2D and 3D' 133 : CALL ERREUR(21) 134 : GOTO 9999 135 : ENDIF 136 : ENDDO 137 : ELSE 138 : LOGCOM=.TRUE. 139 : DO ISOUS=2,NBS,1 140 : IF((MLENTI.LECT(NBS) .EQ. 7) .OR. 141 : & (MLENTI.LECT(NBS) .EQ. 11))THEN 142 : WRITE(IOIMP,*) 'Euler model' 143 : WRITE(IOIMP,*) 'No mixing between 2D and 3D' 143 : CALL ERREUR(21) 144 : GOTO 9999 145 : ENDIF 146 : ENDDO 147 : ENDIF 148 : ENDIF 149 : C 150 : C**** 'MAILLAGE' 151 : C 'CENTRE' 152 : C 153 : C Initialisation 154 : C 155 : SEGINI, MELMAI=MELEME 156 : NBREF=0 157 : IF(NBS .EQ. 1)THEN 158 : ISOUS=1 159 : IF(MLENTI.LECT(ISOUS) .EQ. 3)THEN 160 : C SEG3 -> SEG2 161 : MELMAI.ITYPEL=2 162 : NBNN=2 163 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 7)THEN 164 : C TRI7 -> TRI3 165 : MELMAI.ITYPEL=4 166 : NBNN=3 167 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 11)THEN 168 : C QUA9 -> QUA4 169 : MELMAI.ITYPEL=8 170 : NBNN=4 171 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 35)THEN 172 : C TE15 -> TET4 173 : MELMAI.ITYPEL=23 174 : NBNN=4 175 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 36)THEN 176 : C PY19 -> PYR5 177 : MELMAI.ITYPEL=25 178 : NBNN=5 179 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 34)THEN 180 : C PR21 -> PRI6 181 : MELMAI.ITYPEL=16 182 : NBNN=6 183 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 33)THEN 184 : C CU27 -> CUB8 185 : MELMAI.ITYPEL=14 186 : NBNN=8 187 : ENDIF 188 : NBELEM=MLENT1.LECT(ISOUS) 189 : NBSOUS=0 190 : SEGADJ MELMAI 191 : NBE0=NBELEM 192 : ELSE 193 : NBE0=0 194 : DO ISOUS=1,NBS,1 195 : NBELEM=MLENT1.LECT(ISOUS) 196 : IF(MLENTI.LECT(ISOUS) .EQ. 3)THEN 197 : WRITE(IOIMP,*) 'Subroutine kdom10.eso' 197 : CALL ERREUR(5) 198 : GOTO 9999 199 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 7)THEN 200 : C TRI7 -> TRI3 201 : NBNN=3 202 : MELMAI.ITYPEL=4 203 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 11)THEN 204 : C QUA9 -> QUA4 205 : MELMAI.ITYPEL=8 206 : NBNN=4 207 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 35)THEN 208 : C TE15 -> TET4 209 : MELMAI.ITYPEL=23 210 : NBNN=4 211 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 36)THEN 212 : C PY19 -> PYR5 213 : MELMAI.ITYPEL=25 214 : NBNN=5 215 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 34)THEN 216 : C PR21 -> PRI6 217 : MELMAI.ITYPEL=16 218 : NBNN=6 219 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 33)THEN 220 : C CU27 -> CUB8 221 : MELMAI.ITYPEL=14 222 : NBNN=8 223 : ENDIF 224 : NBSOUS=0 225 : SEGINI IPT2 226 : MELMAI.LISOUS(ISOUS)=IPT2 227 : IPT2.ITYPEL=MELMAI.ITYPEL 228 : MELMAI.ITYPEL=0 229 : NBE0=NBE0+NBELEM 230 : IPT1=MELEME.LISOUS(ISOUS) 231 : ENDDO 232 : ENDIF 233 : C 234 : NBELEM=NBE0 235 : NBNN=1 236 : NBSOUS=0 237 : NBREF=0 238 : SEGINI MELCEN 239 : MELCEN.ITYPEL=1 240 : NBE0=0 241 : C 242 : C**** Filling 243 : C 244 : DO ISOUS=1,NBS,1 245 : IF(NBS.EQ.1)THEN 246 : IPT1=MELEME 247 : IPT2=MELMAI 248 : ELSE 249 : IPT1=MELEME.LISOUS(ISOUS) 250 : IPT2=MELMAI.LISOUS(ISOUS) 251 : ENDIF 252 : NBELEM=MLENT1.LECT(ISOUS) 253 : IF(MLENTI.LECT(ISOUS) .EQ. 3)THEN 254 : C SEG 255 : DO IELEM=1,NBELEM,1 256 : NBE0=NBE0+1 257 : IPT2.NUM(1,IELEM)=IPT1.NUM(1,IELEM) 258 : IPT2.NUM(2,IELEM)=IPT1.NUM(3,IELEM) 259 : IPT2.ICOLOR(IELEM)=IPT1.ICOLOR(IELEM) 260 : MELCEN.NUM(1,NBE0)=IPT1.NUM(2,IELEM) 261 : MELCEN.ICOLOR(NBE0)=IPT1.ICOLOR(IELEM) 262 : ENDDO 263 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 7)THEN 264 : C TRI 265 : DO IELEM=1,NBELEM,1 266 : NBE0=NBE0+1 267 : IPT2.NUM(1,IELEM)=IPT1.NUM(1,IELEM) 268 : IPT2.NUM(2,IELEM)=IPT1.NUM(3,IELEM) 269 : IPT2.NUM(3,IELEM)=IPT1.NUM(5,IELEM) 270 : IPT2.ICOLOR(IELEM)=IPT1.ICOLOR(IELEM) 271 : MELCEN.NUM(1,NBE0)=IPT1.NUM(7,IELEM) 272 : MELCEN.ICOLOR(NBE0)=IPT1.ICOLOR(IELEM) 273 : ENDDO 274 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 11)THEN 275 : C QUA 276 : DO IELEM=1,NBELEM,1 277 : NBE0=NBE0+1 278 : IPT2.NUM(1,IELEM)=IPT1.NUM(1,IELEM) 279 : IPT2.NUM(2,IELEM)=IPT1.NUM(3,IELEM) 280 : IPT2.NUM(3,IELEM)=IPT1.NUM(5,IELEM) 281 : IPT2.NUM(4,IELEM)=IPT1.NUM(7,IELEM) 282 : IPT2.ICOLOR(IELEM)=IPT1.ICOLOR(IELEM) 283 : MELCEN.NUM(1,NBE0)=IPT1.NUM(9,IELEM) 284 : MELCEN.ICOLOR(NBE0)=IPT1.ICOLOR(IELEM) 285 : ENDDO 286 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 35)THEN 287 : C TET 288 : DO IELEM=1,NBELEM,1 289 : NBE0=NBE0+1 290 : IPT2.NUM(1,IELEM)=IPT1.NUM(1,IELEM) 291 : IPT2.NUM(2,IELEM)=IPT1.NUM(3,IELEM) 292 : IPT2.NUM(3,IELEM)=IPT1.NUM(5,IELEM) 293 : IPT2.NUM(4,IELEM)=IPT1.NUM(10,IELEM) 294 : IPT2.ICOLOR(IELEM)=IPT1.ICOLOR(IELEM) 295 : MELCEN.NUM(1,NBE0)=IPT1.NUM(15,IELEM) 296 : MELCEN.ICOLOR(NBE0)=IPT1.ICOLOR(IELEM) 297 : ENDDO 298 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 36)THEN 299 : C PYR 300 : DO IELEM=1,NBELEM,1 301 : NBE0=NBE0+1 302 : IPT2.NUM(1,IELEM)=IPT1.NUM(1,IELEM) 303 : IPT2.NUM(2,IELEM)=IPT1.NUM(3,IELEM) 304 : IPT2.NUM(3,IELEM)=IPT1.NUM(5,IELEM) 305 : IPT2.NUM(4,IELEM)=IPT1.NUM(7,IELEM) 306 : IPT2.NUM(5,IELEM)=IPT1.NUM(13,IELEM) 307 : IPT2.ICOLOR(IELEM)=IPT1.ICOLOR(IELEM) 308 : MELCEN.NUM(1,NBE0)=IPT1.NUM(19,IELEM) 309 : MELCEN.ICOLOR(NBE0)=IPT1.ICOLOR(IELEM) 310 : ENDDO 311 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 34)THEN 312 : C PRI 313 : DO IELEM=1,NBELEM,1 314 : NBE0=NBE0+1 315 : IPT2.NUM(1,IELEM)=IPT1.NUM(1,IELEM) 316 : IPT2.NUM(2,IELEM)=IPT1.NUM(3,IELEM) 317 : IPT2.NUM(3,IELEM)=IPT1.NUM(5,IELEM) 318 : IPT2.NUM(4,IELEM)=IPT1.NUM(10,IELEM) 319 : IPT2.NUM(5,IELEM)=IPT1.NUM(12,IELEM) 320 : IPT2.NUM(6,IELEM)=IPT1.NUM(14,IELEM) 321 : IPT2.ICOLOR(IELEM)=IPT1.ICOLOR(IELEM) 322 : MELCEN.NUM(1,NBE0)=IPT1.NUM(21,IELEM) 323 : MELCEN.ICOLOR(NBE0)=IPT1.ICOLOR(IELEM) 324 : ENDDO 325 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 33)THEN 326 : C CUB 327 : DO IELEM=1,NBELEM,1 328 : NBE0=NBE0+1 329 : IPT2.NUM(1,IELEM)=IPT1.NUM(1,IELEM) 330 : IPT2.NUM(2,IELEM)=IPT1.NUM(3,IELEM) 331 : IPT2.NUM(3,IELEM)=IPT1.NUM(5,IELEM) 332 : IPT2.NUM(4,IELEM)=IPT1.NUM(7,IELEM) 333 : IPT2.NUM(5,IELEM)=IPT1.NUM(13,IELEM) 334 : IPT2.NUM(6,IELEM)=IPT1.NUM(15,IELEM) 335 : IPT2.NUM(7,IELEM)=IPT1.NUM(17,IELEM) 336 : IPT2.NUM(8,IELEM)=IPT1.NUM(19,IELEM) 337 : IPT2.ICOLOR(IELEM)=IPT1.ICOLOR(IELEM) 338 : MELCEN.NUM(1,NBE0)=IPT1.NUM(27,IELEM) 339 : MELCEN.ICOLOR(NBE0)=IPT1.ICOLOR(IELEM) 340 : ENDDO 341 : ELSE 342 : WRITE(IOIMP,*) 'Subroutine kdom10.eso' 342 : CALL ERREUR(5) 343 : ENDIF 344 : IF(NBS .NE. 1)THEN 345 : SEGDES IPT2 346 : ENDIF 347 : ENDDO 347 : CALL ECMO(MTAB,'MAILLAGE','MAILLAGE',MELMAI) 347 : CALL ECMO(MTAB,'CENTRE','MAILLAGE',MELCEN) 348 : SEGDES MELCEN 349 : SEGDES MELMAI 350 : C 351 : C**** Volume 352 : C 353 : TYPI='CENTRE ' 354 : JGN=4 355 : JGM=1 356 : SEGINI MLMOTS 357 : MLMOTS.MOTS(1)='SCAL' 357 : CALL KRCHP1(TYPI,MELCEN,MCHPOI,MLMOTS) 357 : CALL LICHT(MCHPOI,MPOVAL,TYPI,IGEOM) 358 : C SEGACT MPOVAL 359 : C 360 : SEGSUP MLMOTS 361 : NBE0=0 362 : DO ISOUS=1,NBS,1 363 : IF(NBS.EQ.1) THEN 364 : IPT1=MELEME 365 : ELSE 366 : IPT1=MELEME.LISOUS(ISOUS) 367 : ENDIF 368 : C 369 : NBELEM=MLENT1.LECT(ISOUS) 370 : IF(MLENTI.LECT(ISOUS) .EQ. 3)THEN 371 : C SEG3 372 : DO IELEM=1,NBELEM,1 373 : NBE0=NBE0+1 374 : NN1=IPT1.NUM(1,IELEM) 375 : NN2=IPT1.NUM(2,IELEM) 376 : NN3=IPT1.NUM(3,IELEM) 376 : CALL KDOM5(NN1,NN2,NN3,VOL) 377 : MPOVAL.VPOCHA(NBE0,1)=VOL 378 : ENDDO 379 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 7)THEN 380 : C TRI 381 : DO IELEM=1,NBELEM,1 382 : NBE0=NBE0+1 383 : NN1=IPT1.NUM(1,IELEM) 384 : NN2=IPT1.NUM(3,IELEM) 385 : NN3=IPT1.NUM(5,IELEM) 386 : NN4=IPT1.NUM(7,IELEM) 387 : C Vol TRI7 as a vol of a degenerate QUA9 387 : CALL KDOM6(NN1,NN2,NN3,NN3,NN4,VOL) 388 : MPOVAL.VPOCHA(NBE0,1)=VOL 389 : ENDDO 390 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 11)THEN 391 : C QUA 392 : DO IELEM=1,NBELEM,1 393 : NBE0=NBE0+1 394 : NN1=IPT1.NUM(1,IELEM) 395 : NN2=IPT1.NUM(3,IELEM) 396 : NN3=IPT1.NUM(5,IELEM) 397 : NN4=IPT1.NUM(7,IELEM) 398 : NN5=IPT1.NUM(9,IELEM) 398 : CALL KDOM6(NN1,NN2,NN3,NN4,NN5,VOL) 399 : MPOVAL.VPOCHA(NBE0,1)=VOL 400 : ENDDO 401 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 35)THEN 402 : C TET 403 : DO IELEM=1,NBELEM,1 404 : NBE0=NBE0+1 405 : NN1=IPT1.NUM(1,IELEM) 406 : NN2=IPT1.NUM(3,IELEM) 407 : NN3=IPT1.NUM(5,IELEM) 408 : NN4=IPT1.NUM(10,IELEM) 408 : CALL KDOM4(NN1,NN2,NN3,NN4,CEL,VOL) 409 : MPOVAL.VPOCHA(NBE0,1)=ABS(VOL) 410 : ENDDO 411 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 36)THEN 412 : C PYR 413 : DO IELEM=1,NBELEM,1 414 : NBE0=NBE0+1 415 : NN1=IPT1.NUM(1,IELEM) 416 : NN2=IPT1.NUM(3,IELEM) 417 : NN3=IPT1.NUM(5,IELEM) 418 : NN4=IPT1.NUM(7,IELEM) 419 : NN5=IPT1.NUM(14,IELEM) 420 : NN6=IPT1.NUM(13,IELEM) 420 : CALL KDOM3(NN1,NN2,NN3,NN4,NN5,NN6,CEL,VOL) 421 : MPOVAL.VPOCHA(NBE0,1)=ABS(VOL) 422 : ENDDO 423 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 34)THEN 424 : C PRI 425 : C To compute the volume, we divide it into 2 426 : C TET et 3 PYR 427 : C 428 : DO IELEM=1,NBELEM,1 429 : VOL=0.0D0 430 : NBE0=NBE0+1 431 : NN1=IPT1.NUM(1,IELEM) 432 : NN2=IPT1.NUM(3,IELEM) 433 : NN3=IPT1.NUM(5,IELEM) 434 : NN4=IPT1.NUM(21,IELEM) 434 : CALL KDOM4(NN1,NN2,NN3,NN4,CEL,VOL1) 435 : VOL=VOL+VOL1 436 : NN1=IPT1.NUM(14,IELEM) 437 : NN2=IPT1.NUM(12,IELEM) 438 : NN3=IPT1.NUM(10,IELEM) 439 : NN4=IPT1.NUM(21,IELEM) 439 : CALL KDOM4(NN1,NN2,NN3,NN4,CEL,VOL1) 440 : VOL=VOL+VOL1 441 : NN1=IPT1.NUM(1,IELEM) 442 : NN2=IPT1.NUM(10,IELEM) 443 : NN3=IPT1.NUM(12,IELEM) 444 : NN4=IPT1.NUM(3,IELEM) 445 : NN5=IPT1.NUM(16,IELEM) 446 : NN6=IPT1.NUM(21,IELEM) 446 : CALL KDOM3(NN1,NN2,NN3,NN4,NN5,NN6,CEL,VOL1) 447 : VOL=VOL+VOL1 448 : NN1=IPT1.NUM(3,IELEM) 449 : NN2=IPT1.NUM(12,IELEM) 450 : NN3=IPT1.NUM(14,IELEM) 451 : NN4=IPT1.NUM(5,IELEM) 452 : NN5=IPT1.NUM(17,IELEM) 453 : NN6=IPT1.NUM(21,IELEM) 453 : CALL KDOM3(NN1,NN2,NN3,NN4,NN5,NN6,CEL,VOL1) 454 : VOL=VOL+VOL1 455 : NN1=IPT1.NUM(10,IELEM) 456 : NN2=IPT1.NUM(1,IELEM) 457 : NN3=IPT1.NUM(5,IELEM) 458 : NN4=IPT1.NUM(14,IELEM) 459 : NN5=IPT1.NUM(18,IELEM) 460 : NN6=IPT1.NUM(21,IELEM) 460 : CALL KDOM3(NN1,NN2,NN3,NN4,NN5,NN6,CEL,VOL1) 461 : VOL=VOL+VOL1 462 : MPOVAL.VPOCHA(NBE0,1)=ABS(VOL) 463 : ENDDO 464 : C 465 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 33)THEN 466 : C CUB 467 : C To compute the volume, we divide it into 6 PYR 468 : C 469 : DO IELEM=1,NBELEM,1 470 : VOL=0.0D0 471 : NBE0=NBE0+1 472 : NN1=IPT1.NUM(1,IELEM) 473 : NN2=IPT1.NUM(3,IELEM) 474 : NN3=IPT1.NUM(5,IELEM) 475 : NN4=IPT1.NUM(7,IELEM) 476 : NN5=IPT1.NUM(25,IELEM) 477 : NN6=IPT1.NUM(27,IELEM) 477 : CALL KDOM3(NN1,NN2,NN3,NN4,NN5,NN6,CEL,VOL1) 478 : VOL=VOL+VOL1 479 : NN1=IPT1.NUM(19,IELEM) 480 : NN2=IPT1.NUM(17,IELEM) 481 : NN3=IPT1.NUM(15,IELEM) 482 : NN4=IPT1.NUM(13,IELEM) 483 : NN5=IPT1.NUM(26,IELEM) 484 : NN6=IPT1.NUM(27,IELEM) 484 : CALL KDOM3(NN1,NN2,NN3,NN4,NN5,NN6,CEL,VOL1) 485 : VOL=VOL+VOL1 486 : NN1=IPT1.NUM(1,IELEM) 487 : NN2=IPT1.NUM(13,IELEM) 488 : NN3=IPT1.NUM(15,IELEM) 489 : NN4=IPT1.NUM(3,IELEM) 490 : NN5=IPT1.NUM(21,IELEM) 491 : NN6=IPT1.NUM(27,IELEM) 491 : CALL KDOM3(NN1,NN2,NN3,NN4,NN5,NN6,CEL,VOL1) 492 : VOL=VOL+VOL1 493 : NN1=IPT1.NUM(3,IELEM) 494 : NN2=IPT1.NUM(15,IELEM) 495 : NN3=IPT1.NUM(17,IELEM) 496 : NN4=IPT1.NUM(5,IELEM) 497 : NN5=IPT1.NUM(22,IELEM) 498 : NN6=IPT1.NUM(27,IELEM) 498 : CALL KDOM3(NN1,NN2,NN3,NN4,NN5,NN6,CEL,VOL1) 499 : VOL=VOL+VOL1 500 : NN1=IPT1.NUM(7,IELEM) 501 : NN2=IPT1.NUM(5,IELEM) 502 : NN3=IPT1.NUM(17,IELEM) 503 : NN4=IPT1.NUM(19,IELEM) 504 : NN5=IPT1.NUM(23,IELEM) 505 : NN6=IPT1.NUM(27,IELEM) 505 : CALL KDOM3(NN1,NN2,NN3,NN4,NN5,NN6,CEL,VOL1) 506 : VOL=VOL+VOL1 507 : NN1=IPT1.NUM(19,IELEM) 508 : NN2=IPT1.NUM(13,IELEM) 509 : NN3=IPT1.NUM(1,IELEM) 510 : NN4=IPT1.NUM(7,IELEM) 511 : NN5=IPT1.NUM(24,IELEM) 512 : NN6=IPT1.NUM(27,IELEM) 512 : CALL KDOM3(NN1,NN2,NN3,NN4,NN5,NN6,CEL,VOL1) 513 : VOL=VOL+VOL1 514 : MPOVAL.VPOCHA(NBE0,1)=ABS(VOL) 515 : ENDDO 516 : ELSE 517 : WRITE(IOIMP,*) 'Subroutine kdom10.eso' 517 : CALL ERREUR(5) 518 : ENDIF 519 : ENDDO 520 : C 521 : SEGDES MPOVAL 521 : CALL ECMO(MTAB,'XXVOLUM','CHPOINT',MCHPOI) 522 : C 523 : C 524 : C**** For the moment we have: 525 : C 'QUAF' 526 : C 'MAILLAGE' 527 : C 'CENTRE' 528 : C 'XXVOLUM' 529 : C 530 : C In the complete case (2D mesh in 2D and 3D mesh in 3D): 531 : C 532 : IF(LOGCOM)THEN 533 : C 534 : C******* We create ELTFA 535 : C 536 : SEGINI, MELTFA=MELEME 537 : NBREF=0 538 : IF(NBS .EQ. 1)THEN 539 : ISOUS=1 540 : IF(MLENTI.LECT(ISOUS) .EQ. 7)THEN 541 : C TRI3 542 : NBNN=3 543 : MELTFA.ITYPEL=4 544 : C ELTFA TRI3 545 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 11)THEN 546 : C QUA4 547 : NBNN=4 548 : MELTFA.ITYPEL=8 549 : C ELTFA QUA4 550 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 35)THEN 551 : C TET4 552 : NBNN=4 553 : MELTFA.ITYPEL=23 554 : C ELTFA TET4 555 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 36)THEN 556 : C PYR5 557 : NBNN=5 558 : MELTFA.ITYPEL=9 559 : C ELTFA QUA5 560 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 34)THEN 561 : C PRI6 562 : NBNN=5 563 : MELTFA.ITYPEL=25 564 : C ELTFA PYR5 565 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 33)THEN 566 : C CUB8 567 : NBNN=6 568 : MELTFA.ITYPEL=16 569 : C ELTFA PRI6 570 : ELSE 571 : WRITE(IOIMP,*) 'Subroutine kdom10.eso' 571 : CALL ERREUR(5) 572 : GOTO 9999 573 : ENDIF 574 : NBELEM=MLENT1.LECT(ISOUS) 575 : NBSOUS=0 576 : SEGADJ MELTFA 577 : ELSE 578 : DO ISOUS=1,NBS,1 579 : NBELEM=MLENT1.LECT(ISOUS) 580 : IF(MLENTI.LECT(ISOUS) .EQ. 7)THEN 581 : C TRI3 582 : NBNN=3 583 : MELTFA.ITYPEL=4 584 : C ELTFA TRI3 585 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 11)THEN 586 : C QUA4 587 : NBNN=4 588 : MELTFA.ITYPEL=8 589 : C ELTFA QUA4 590 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 35)THEN 591 : C TET4 592 : NBNN=4 593 : MELTFA.ITYPEL=23 594 : C ELTFA TET4 595 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 36)THEN 596 : C PYR5 597 : NBNN=5 598 : MELTFA.ITYPEL=9 599 : C ELTFA QUA5 600 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 34)THEN 601 : C PRI6 602 : NBNN=5 603 : MELTFA.ITYPEL=25 604 : C ELTFA PYR5 605 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 33)THEN 606 : C CUB8 607 : NBNN=6 608 : MELTFA.ITYPEL=16 609 : C ELTFA PRI6 610 : ELSE 611 : WRITE(IOIMP,*) 'Subroutine kdom10.eso' 611 : CALL ERREUR(5) 612 : GOTO 9999 613 : ENDIF 614 : NBSOUS=0 615 : SEGINI IPT2 616 : MELTFA.LISOUS(ISOUS)=IPT2 617 : C 618 : IPT2.ITYPEL=MELTFA.ITYPEL 619 : MELTFA.ITYPEL=0 620 : ENDDO 621 : ENDIF 622 : C 623 : C******* We fill ELTFA 624 : C We also count: 625 : C NFACT = number of triangular faces 626 : C NFAC = number of non-triangular faces 627 : C NSOM = number of SOMMET 628 : C 629 : NTP=MCOORD.XCOOR(/1)/(IDIM+1) 630 : JG=NTP 631 : NFAC=0 632 : NFACT=0 633 : NSOM=0 634 : 635 : LAST=-1 636 : SEGINI MLRES 637 : LASTT=-1 638 : SEGINI MLREST 639 : LASTS=-1 640 : SEGINI MLRESS 641 : C LAST+MLRES = chaining list to find the (non-triangular faces) 642 : DO ISOUS=1,NBS,1 643 : IF(NBS.EQ.1) THEN 644 : IPT1=MELEME 645 : IPT2=MELTFA 646 : ELSE 647 : IPT1=MELEME.LISOUS(ISOUS) 648 : IPT2=MELTFA.LISOUS(ISOUS) 649 : ENDIF 650 : C 651 : NBELEM=MLENT1.LECT(ISOUS) 652 : IF(MLENTI.LECT(ISOUS) .EQ. 7)THEN 653 : C TRI (2D) 654 : LISFAC(1)=2 655 : LISFAC(2)=4 656 : LISFAC(3)=6 657 : LISSOM(1)=1 658 : LISSOM(2)=3 659 : LISSOM(3)=5 660 : NNS=3 661 : NNOEU=3 662 : NNOEUT=0 663 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 11)THEN 664 : C QUA (2D) 665 : LISFAC(1)=2 666 : LISFAC(2)=4 667 : LISFAC(3)=6 668 : LISFAC(4)=8 669 : LISSOM(1)=1 670 : LISSOM(2)=3 671 : LISSOM(3)=5 672 : LISSOM(4)=7 673 : NNS=4 674 : NNOEU=4 675 : NNOEUT=0 676 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 35)THEN 677 : C TET 678 : LISFAT(1)=11 679 : LISFAT(2)=12 680 : LISFAT(3)=13 681 : LISFAT(4)=14 682 : LISSOM(1)=1 683 : LISSOM(2)=3 684 : LISSOM(3)=5 685 : LISSOM(4)=10 686 : NNS=4 687 : NNOEU=0 688 : NNOEUT=4 689 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 36)THEN 690 : C PYR 691 : LISFAC(1)=14 692 : LISFAT(1)=15 693 : LISFAT(2)=16 694 : LISFAT(3)=17 695 : LISFAT(4)=18 696 : LISSOM(1)=1 697 : LISSOM(2)=3 698 : LISSOM(3)=5 699 : LISSOM(4)=7 700 : LISSOM(5)=13 701 : NNS=5 702 : NNOEU=1 703 : NNOEUT=4 704 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 34)THEN 705 : C PRI 706 : LISFAC(1)=16 707 : LISFAC(2)=17 708 : LISFAC(3)=18 709 : LISFAT(1)=19 710 : LISFAT(2)=20 711 : LISSOM(1)=1 712 : LISSOM(2)=3 713 : LISSOM(3)=5 714 : LISSOM(4)=10 715 : LISSOM(5)=12 716 : LISSOM(6)=14 717 : NNS=6 718 : NNOEU=3 719 : NNOEUT=2 720 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 33)THEN 721 : C CUB 722 : LISFAC(1)=25 723 : LISFAC(2)=26 724 : LISFAC(3)=21 725 : LISFAC(4)=22 726 : LISFAC(5)=23 727 : LISFAC(6)=24 728 : LISSOM(1)=1 729 : LISSOM(2)=3 730 : LISSOM(3)=5 731 : LISSOM(4)=7 732 : LISSOM(5)=13 733 : LISSOM(6)=15 734 : LISSOM(7)=17 735 : LISSOM(8)=19 736 : NNS=8 737 : NNOEU=6 738 : NNOEUT=0 739 : ELSE 740 : WRITE(IOIMP,*) 'Subroutine kdom10.eso' 740 : CALL ERREUR(5) 741 : ENDIF 742 : C 743 : DO IELEM=1,NBELEM,1 744 : IFAC=0 745 : DO INOE=1,NNOEU,1 746 : NN1=IPT1.NUM(LISFAC(INOE),IELEM) 747 : IFAC=IFAC+1 748 : IPT2.NUM(IFAC,IELEM)=NN1 749 : IF(MLRES.LECT(NN1) .EQ. 0)THEN 750 : NFAC=NFAC+1 751 : MLRES.LECT(NN1)=LAST 752 : LAST=NN1 753 : ENDIF 754 : ENDDO 755 : DO INOE=1,NNOEUT,1 756 : NN1=IPT1.NUM(LISFAT(INOE),IELEM) 757 : IFAC=IFAC+1 758 : IPT2.NUM(IFAC,IELEM)=NN1 759 : IF(MLREST.LECT(NN1) .EQ. 0)THEN 760 : NFACT=NFACT+1 761 : MLREST.LECT(NN1)=LASTT 762 : LASTT=NN1 763 : ENDIF 764 : ENDDO 765 : DO INOE=1,NNS,1 766 : NN1=IPT1.NUM(LISSOM(INOE),IELEM) 767 : IF(MLRESS.LECT(NN1) .EQ. 0)THEN 768 : NSOM=NSOM+1 769 : MLRESS.LECT(NN1)=LASTS 770 : LASTS=NN1 771 : ENDIF 772 : ENDDO 773 : C 774 : C************* Cas particulier: PRISME 775 : C 776 : IF(MLENTI.LECT(ISOUS) .EQ. 34)THEN 777 : NN1=IPT2.NUM(1,IELEM) 778 : NN2=IPT2.NUM(2,IELEM) 779 : NN3=IPT2.NUM(3,IELEM) 780 : DO INOE=1,2 781 : IPT2.NUM(INOE,IELEM)=IPT2.NUM(INOE+3,IELEM) 782 : ENDDO 783 : IPT2.NUM(3,IELEM)=NN1 784 : IPT2.NUM(4,IELEM)=NN2 785 : IPT2.NUM(5,IELEM)=NN3 786 : ENDIF 787 : ENDDO 788 : C 789 : IF(NBS.NE.1) THEN 790 : SEGDES IPT2 791 : ENDIF 792 : C 793 : ENDDO 794 : SEGDES MELTFA 794 : CALL ECMO(MTAB,'ELTFA','MAILLAGE',MELTFA) 795 : 796 : CC 797 : CC******* Test 798 : CC 799 : C write(*,*) 'Triangular faces ', nfact 800 : C 10 if(lastt .ne. -1)then 801 : C write(*,*) lastt 802 : C lastt=mlrest.lect(lastt) 803 : C endif 804 : C if(lastt .ne. -1) goto 10 805 : C write(*,*) 'Non triangular faces ', nfac 806 : C 20 if(last .ne. -1)then 807 : C write(*,*) last 808 : C last=mlres.lect(last) 809 : C endif 810 : C if(last .ne. -1) goto 20 811 : C write(*,*) 'Sommets ', nsom 812 : C 30 if(lasts .ne. -1)then 813 : C write(*,*) lasts 814 : C lasts=mlress.lect(lasts) 815 : C endif 816 : C if(lasts .ne. -1) goto 30 817 : C 818 : CC 819 : CC******* Fin test 820 : CC 821 : C 822 : C******** Creation of SOMMET 823 : C 824 : NBELEM=NSOM 825 : NBNN=1 826 : NBSOUS=0 827 : NBREF=0 828 : SEGINI MELSOM 829 : MELSOM.ITYPEL=1 830 : DO IELEM=1,NSOM,1 831 : MELSOM.NUM(1,IELEM)=LASTS 832 : LASTS=MLRESS.LECT(LASTS) 833 : ENDDO 834 : IF(LASTS .NE. -1)THEN 835 : WRITE(IOIMP,*) 'Subroutine kdom10.eso' 835 : CALL ERREUR(5) 836 : ENDIF 837 : SEGDES MELSOM 837 : CALL ECMO(MTAB,'SOMMET','MAILLAGE',MELSOM) 838 : SEGSUP MLRESS 839 : C 840 : C******** Creation of FACE (in 3D, triangle first) 841 : C 842 : NBELEM=NFAC+NFACT 843 : NBNN=1 844 : NBSOUS=0 845 : NBREF=0 846 : SEGINI MELFAC 847 : MELFAC.ITYPEL=1 848 : DO IELEM=1,NFACT,1 849 : MELFAC.NUM(1,IELEM)=LASTT 850 : LASTT=MLREST.LECT(LASTT) 851 : ENDDO 852 : IF(LASTT .NE. -1)THEN 853 : WRITE(IOIMP,*) 'Subroutine kdom10.eso' 853 : CALL ERREUR(5) 854 : ENDIF 855 : DO IELEM=1,NFAC,1 856 : MELFAC.NUM(1,NFACT+IELEM)=LAST 857 : LAST=MLRES.LECT(LAST) 858 : ENDDO 859 : IF(LAST .NE. -1)THEN 860 : WRITE(IOIMP,*) 'Subroutine kdom10.eso' 860 : CALL ERREUR(5) 861 : ENDIF 862 : SEGDES MELFAC 863 : SEGSUP MLRES 864 : SEGSUP MLREST 864 : CALL ECMO(MTAB,'FACE','MAILLAGE',MELFAC) 865 : C 866 : C******* Creation of FACEL and FACEP 867 : C 867 : CALL KRIPAD(MELFAC,MLEFAC) 868 : C SEGINI MLEFAC 869 : JG=NFAC+NFACT 870 : SEGINI MLETOF 871 : C 872 : C MLETOF.LECT(I1) = how many times has the i-th face of MELFAP 873 : C already been touched? 874 : C 875 : NBELEM=NFAC+NFACT 876 : NBNN=3 877 : NBSOUS=0 878 : NBREF=0 879 : SEGINI MELFAL 880 : MELFAL.ITYPEL=3 881 : C 882 : IF(IDIM.EQ.2)THEN 883 : C FACEP is a SEG3 884 : NBELEM=NFAC 885 : NBNN=3 886 : NBSOUS=0 887 : NBREF=0 888 : SEGINI MELFAP 889 : MELFAP.ITYPEL=3 890 : IPTQ=MELFAP 891 : ELSEIF(NFAC.EQ.0)THEN 892 : C In 3D we have triangles only 893 : C TRI4 894 : NBELEM=NFACT 895 : NBSOUS=0 896 : NBREF=0 897 : NBNN=4 898 : SEGINI MELFAP 899 : IPTT=MELFAP 900 : MELFAP.ITYPEL=5 901 : ELSEIF(NFACT.EQ.0)THEN 902 : C In 3D we have quadrangles only 903 : C QUA5 904 : NBELEM=NFAC 905 : NBSOUS=0 906 : NBREF=0 907 : NBNN=5 908 : SEGINI MELFAP 909 : IPTQ=MELFAP 910 : MELFAP.ITYPEL=9 911 : ELSE 912 : C TRI4 5 913 : C QUA5 9 914 : NBELEM=0 915 : NBNN=0 916 : NBSOUS=2 917 : NBREF=0 918 : SEGINI MELFAP 919 : MELFAP.ITYPEL=0 920 : C 921 : NBELEM=NFACT 922 : NBSOUS=0 923 : NBREF=0 924 : NBNN=4 925 : SEGINI IPTT 926 : MELFAP.LISOUS(1)=IPTT 927 : IPTT.ITYPEL=5 928 : C 929 : NBELEM=NFAC 930 : NBSOUS=0 931 : NBREF=0 932 : NBNN=5 933 : SEGINI IPTQ 934 : MELFAP.LISOUS(2)=IPTQ 935 : IPTQ.ITYPEL=9 936 : ENDIF 937 : C NBSFP = NBSOUS for FACEP 938 : C 939 : NBSFP=NBSOUS 940 : C 941 : DO ISOUS=1,NBS,1 942 : C 943 : C********** Loop on the elementary mesh of the QUAF 944 : C 944 : C We recall that 945 : C MLENTI.LECT(ISOUS) contains the type of support of 946 : C MELEME.LISOUS 947 : C MLENT1.LECT(ISOUS) contains the number of elements of 948 : C MELEME.LISOUS 949 : C 950 : IF(NBS.EQ.1) THEN 951 : IPT1=MELEME 952 : ELSE 953 : IPT1=MELEME.LISOUS(ISOUS) 954 : ENDIF 955 : C 956 : IF(MLENTI.LECT(ISOUS) .EQ. 7)THEN 957 : C TRI (2D) 958 : LIFAC(1,1)=2 959 : LIFAC(2,1)=1 960 : LIFAC(3,1)=3 961 : LIFAC(1,2)=4 962 : LIFAC(2,2)=3 963 : LIFAC(3,2)=5 964 : LIFAC(1,3)=6 965 : LIFAC(2,3)=5 966 : LIFAC(3,3)=1 967 : C Here we put the center point in LISSOM 968 : LISSOM(1)=7 969 : C Number of non-triangular and triangular faces 970 : NNOEU=3 971 : NNOEUT=0 972 : C 973 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 11)THEN 974 : C QUA (2D) 975 : LIFAC(1,1)=2 976 : LIFAC(2,1)=1 977 : LIFAC(3,1)=3 978 : LIFAC(1,2)=4 979 : LIFAC(2,2)=3 980 : LIFAC(3,2)=5 981 : LIFAC(1,3)=6 982 : LIFAC(2,3)=5 983 : LIFAC(3,3)=7 984 : LIFAC(1,4)=8 985 : LIFAC(2,4)=7 986 : LIFAC(3,4)=1 987 : LISSOM(1)=9 988 : NNOEU=4 989 : NNOEUT=0 990 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 35)THEN 991 : C TET 992 : LIFACT(1,1)=11 993 : LIFACT(2,1)=1 994 : LIFACT(3,1)=3 995 : LIFACT(4,1)=5 996 : LIFACT(1,2)=12 997 : LIFACT(2,2)=1 998 : LIFACT(3,2)=10 999 : LIFACT(4,2)=3 1000 : LIFACT(1,3)=13 1001 : LIFACT(2,3)=3 1002 : LIFACT(3,3)=10 1003 : LIFACT(4,3)=5 1004 : LIFACT(1,4)=14 1005 : LIFACT(2,4)=10 1006 : LIFACT(3,4)=1 1007 : LIFACT(4,4)=5 1008 : LISSOM(1)=15 1009 : NNOEU=0 1010 : NNOEUT=4 1011 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 36)THEN 1012 : C PYR 1013 : LIFAC(1,1)=14 1014 : LIFAC(2,1)=1 1015 : LIFAC(3,1)=3 1016 : LIFAC(4,1)=5 1017 : LIFAC(5,1)=7 1018 : LIFACT(1,1)=15 1019 : LIFACT(2,1)=1 1020 : LIFACT(3,1)=13 1021 : LIFACT(4,1)=3 1022 : LIFACT(1,2)=16 1023 : LIFACT(2,2)=3 1024 : LIFACT(3,2)=13 1025 : LIFACT(4,2)=5 1026 : LIFACT(1,3)=17 1027 : LIFACT(2,3)=5 1028 : LIFACT(3,3)=13 1029 : LIFACT(4,3)=7 1030 : LIFACT(1,4)=18 1031 : LIFACT(2,4)=13 1032 : LIFACT(3,4)=1 1033 : LIFACT(4,4)=7 1034 : LISSOM(1)=19 1035 : NNOEU=1 1036 : NNOEUT=4 1037 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 34)THEN 1038 : C PRI 1039 : LIFAC(1,1)=16 1040 : LIFAC(2,1)=1 1041 : LIFAC(3,1)=10 1042 : LIFAC(4,1)=12 1043 : LIFAC(5,1)=3 1044 : LIFAC(1,2)=17 1045 : LIFAC(2,2)=3 1046 : LIFAC(3,2)=12 1047 : LIFAC(4,2)=14 1048 : LIFAC(5,2)=5 1049 : LIFAC(1,3)=18 1050 : LIFAC(2,3)=1 1051 : LIFAC(3,3)=5 1052 : LIFAC(4,3)=14 1053 : LIFAC(5,3)=10 1054 : LIFACT(1,1)=19 1055 : LIFACT(2,1)=1 1056 : LIFACT(3,1)=3 1057 : LIFACT(4,1)=5 1058 : LIFACT(1,2)=20 1059 : LIFACT(2,2)=10 1060 : LIFACT(3,2)=14 1061 : LIFACT(4,2)=12 1062 : LISSOM(1)=21 1063 : NNOEU=3 1064 : NNOEUT=2 1065 : ELSEIF(MLENTI.LECT(ISOUS) .EQ. 33)THEN 1066 : C CUB 1067 : LIFAC(1,1)=21 1068 : LIFAC(2,1)=1 1069 : LIFAC(3,1)=13 1070 : LIFAC(4,1)=15 1071 : LIFAC(5,1)=3 1072 : LIFAC(1,2)=22 1073 : LIFAC(2,2)=3 1074 : LIFAC(3,2)=15 1075 : LIFAC(4,2)=17 1076 : LIFAC(5,2)=5 1077 : LIFAC(1,3)=23 1078 : LIFAC(2,3)=5 1079 : LIFAC(3,3)=17 1080 : LIFAC(4,3)=19 1081 : LIFAC(5,3)=7 1082 : LIFAC(1,4)=24 1083 : LIFAC(2,4)=1 1084 : LIFAC(3,4)=7 1085 : LIFAC(4,4)=19 1086 : LIFAC(5,4)=13 1087 : LIFAC(1,5)=25 1088 : LIFAC(2,5)=1 1089 : LIFAC(3,5)=3 1090 : LIFAC(4,5)=5 1091 : LIFAC(5,5)=7 1092 : LIFAC(1,6)=26 1093 : LIFAC(2,6)=13 1094 : LIFAC(3,6)=19 1095 : LIFAC(4,6)=17 1096 : LIFAC(5,6)=15 1097 : LISSOM(1)=27 1098 : NNOEU=6 1099 : NNOEUT=0 1100 : ELSE 1101 : WRITE(IOIMP,*) 'Subroutine kdom10.eso' 1101 : CALL ERREUR(5) 1102 : ENDIF 1103 : C 1104 : NBELEM=MLENT1.LECT(ISOUS) 1105 : DO IELEM=1,NBELEM,1 1106 : C NNOEU = number of quagrangular elements 1107 : DO INOE=1,NNOEU,1 1108 : C NN1 is the global number of the face 1109 : C NN2 is the local number of the face in the MELEME 1110 : C 'FACE' 1111 : C The MELEME 'FACE' contains : -triangular faces 1112 : C (MELEME IPTT) 1113 : C -non-triangular faces 1114 : C (MELEME IPTQ) 1115 : C Then NN3 = position of NN1 in IPTQ = NN2 - NFACT 1116 : C where NFACT is the total number of triangular faces 1117 : C 1118 : NN1=IPT1.NUM(LIFAC(1,INOE),IELEM) 1119 : NN2=MLEFAC.LECT(NN1) 1120 : NN3=NN2-NFACT 1121 : IF(MLETOF.LECT(NN2).EQ.0)THEN 1122 : C 1123 : C MLETOF.LECT(NN2) = how many times the face NN2 has 1124 : C been touched? 1125 : C 1126 : MLETOF.LECT(NN2)=1 1127 : MELFAL.NUM(2,NN2)=NN1 1128 : MELFAL.NUM(1,NN2)=IPT1.NUM(LISSOM(1),IELEM) 1129 : MELFAL.NUM(3,NN2)=IPT1.NUM(LISSOM(1),IELEM) 1130 : IF(IDIM .EQ.2)THEN 1131 : IPTQ.NUM(1,NN3)=IPT1.NUM(LIFAC(2,INOE),IELEM) 1132 : IPTQ.NUM(2,NN3)=IPT1.NUM(LIFAC(3,INOE),IELEM) 1133 : IPTQ.NUM(3,NN3)=NN1 1134 : ELSE 1135 : IPTQ.NUM(1,NN3)=IPT1.NUM(LIFAC(2,INOE),IELEM) 1136 : IPTQ.NUM(2,NN3)=IPT1.NUM(LIFAC(3,INOE),IELEM) 1137 : IPTQ.NUM(3,NN3)=IPT1.NUM(LIFAC(4,INOE),IELEM) 1138 : IPTQ.NUM(4,NN3)=IPT1.NUM(LIFAC(5,INOE),IELEM) 1139 : IPTQ.NUM(5,NN3)=NN1 1140 : ENDIF 1141 : ELSEIF(MLETOF.LECT(NN2).EQ.1)THEN 1142 : MLETOF.LECT(NN2)=2 1143 : MELFAL.NUM(3,NN2)=IPT1.NUM(LISSOM(1),IELEM) 1144 : ELSE 1145 : WRITE(IOIMP,*) 'Subroutine kdom10.eso' 1145 : CALL ERREUR(5) 1146 : ENDIF 1147 : ENDDO 1148 : DO INOE=1,NNOEUT,1 1149 : C NN1 is the global number of the face 1150 : C NN2 is the local number of the face in the MELEME 1151 : C 'FACE' and is also the position of NN1 in IPTT 1152 : C 1153 : NN1=IPT1.NUM(LIFACT(1,INOE),IELEM) 1154 : NN2=MLEFAC.LECT(NN1) 1155 : IF(MLETOF.LECT(NN2).EQ.0)THEN 1156 : MLETOF.LECT(NN2)=1 1157 : MELFAL.NUM(2,NN2)=NN1 1158 : MELFAL.NUM(1,NN2)=IPT1.NUM(LISSOM(1),IELEM) 1159 : MELFAL.NUM(3,NN2)=IPT1.NUM(LISSOM(1),IELEM) 1160 : IPTT.NUM(1,NN2)=IPT1.NUM(LIFACT(2,INOE),IELEM) 1161 : IPTT.NUM(2,NN2)=IPT1.NUM(LIFACT(3,INOE),IELEM) 1162 : IPTT.NUM(3,NN2)=IPT1.NUM(LIFACT(4,INOE),IELEM) 1163 : IPTT.NUM(4,NN2)=NN1 1164 : ELSEIF(MLETOF.LECT(NN2).EQ.1)THEN 1165 : MLETOF.LECT(NN2)=2 1166 : MELFAL.NUM(3,NN2)=IPT1.NUM(LISSOM(1),IELEM) 1167 : ELSE 1168 : WRITE(IOIMP,*) 'Subroutine kdom10.eso' 1168 : CALL ERREUR(5) 1169 : ENDIF 1170 : ENDDO 1171 : ENDDO 1172 : ENDDO 1173 : IF(NBSFP .NE. 0)THEN 1174 : SEGDES IPTQ 1175 : SEGDES IPTT 1176 : ENDIF 1177 : SEGDES MELFAL 1178 : SEGDES MELFAP 1179 : C SEGDES MELFAP 1180 : C SEGDES MELFAP 1181 : C SEGDES MELFAP 1182 : SEGSUP MLETOF 1183 : SEGSUP MLEFAC 1183 : CALL ECMO(MTAB,'FACEL','MAILLAGE',MELFAL) 1183 : CALL ECMO(MTAB,'FACEP','MAILLAGE',MELFAP) 1184 : C 1185 : C******** Creation of 'XXSURFAC', 'XXNORMAF', 'MATROT' 1186 : C 1186 : CALL KDOM11(MELFAC,MELFAL,MELFAP,MCHPSU,MCHPNO,MCHPMR) 1187 : IF(IERR.NE.0)GOTO 9999 1187 : CALL ECMO(MTAB,'XXSURFAC','CHPOINT',MCHPSU) 1187 : CALL ECMO(MTAB,'XXNORMAF','CHPOINT',MCHPNO) 1187 : CALL ECMO(MTAB,'MATROT','CHPOINT',MCHPMR) 1188 : C 1189 : C******** Creation of 'XXDIEMIN' 1190 : C 1190 : CALL KDOM12(MELTFA,MELCEN,MELFAC,MCHPNO,MCHDIA) 1191 : IF(IERR.NE.0)GOTO 9999 1191 : CALL ECMO(MTAB,'XXDIEMIN','CHPOINT',MCHDIA) 1192 : C 1193 : ENDIF 1194 : C 1195 : C 1196 : C***** Fin qui 1197 : C 1198 : IF(NBS .NE. 1)THEN 1199 : DO ISOUS=1,NBS,1 1200 : IPT1=MELEME.LISOUS(ISOUS) 1201 : SEGDES IPT1 1202 : ENDDO 1203 : ENDIF 1204 : C 1205 : SEGSUP MLENTI 1206 : SEGSUP MLENT1 1207 : SEGDES MELEME 1208 : C 1209 : 9999 RETURN 1210 : C 1211 : END 1212 : 1213 :
© Cast3M 2003 - All rights reserved.
Disclaimer