1 : C KDOM2 SOURCE CHAT 05/01/13 00:53:53 5004 2 : SUBROUTINE KDOM2(MELEMQ) 3 : C 4 : C************************************************************************ 5 : C 6 : C PROJET : CASTEM 2000 7 : C 8 : C NOM : KDOM2 9 : C 9 : C DESCRIPTION : Subroutine called by KDOM1 10 : C In the case of EULER model, we change the 11 : C position of the QUAF points 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 MELEMQ, NBSOUS, ISOUS, NBEL, IEL, NN1, NN2, NN3, I1 35 : & ,NN4,NN5,NN6,NN7,NN8,NN9,NN10,NN11,NN12,NN13,NN14,NN15 36 : & ,NN16,NN17,NN18,NN19,NN20,NN21,NN22,NN23,NN24,NN25,NN26 37 : & ,NN27 38 : REAL*8 XCEN(3),VOL,X1(3),X2(3),X3(3),X4(3),X5(3),X6(3) 39 : & ,VOL1,VOL2,VOL3,VOL4,VOL5,VOL6 40 : LOGICAL LOSEG3, LOTRI7, LOQUA9, LOTE15, LOPY19, LOPR21, LOCU27 41 : C 42 : -INC CCOPTIO 43 : -INC SMELEME 44 : -INC SMCOORD 45 : C 46 : C Elements allowed are SEG3, TRI7, QUA9, TE15, 'PY19', 'PR21', 'CU27' 47 : C ITYPEL 3 7 11 35 36 34 33 48 : C They can be referred just once 49 : C 50 : LOSEG3=.FALSE. 51 : LOTRI7=.FALSE. 52 : LOQUA9=.FALSE. 53 : LOTE15=.FALSE. 54 : LOPY19=.FALSE. 55 : LOPR21=.FALSE. 56 : LOCU27=.FALSE. 57 : C 58 : MELEME=MELEMQ 59 : SEGACT MELEME 60 : NBSOUS=MELEME.LISOUS(/1) 61 : IF(NBSOUS .EQ. 0) NBSOUS=1 62 : DO ISOUS=1,NBSOUS,1 63 : IF(NBSOUS .NE. 1)THEN 64 : IPT1=MELEME.LISOUS(ISOUS) 65 : SEGACT IPT1 66 : ELSE 67 : IPT1=MELEME 68 : ENDIF 69 : NBEL=IPT1.NUM(/2) 70 : C 71 : IF(IPT1.ITYPEL .EQ. 3)THEN 72 : IF(LOSEG3)THEN 73 : C Elt already referred 74 : WRITE(IOIMP,*) 'Subroutine kdom2' 75 : WRITE(IOIMP,*) 'Mesh type not recognized' 75 : CALL ERREUR(5) 76 : ENDIF 77 : LOSEG3=.TRUE. 78 : C 79 : C********** SEG3 80 : C 81 : DO IEL=1,NBEL,1 82 : NN1=IPT1.NUM(1,IEL) 83 : NN2=IPT1.NUM(2,IEL) 84 : NN3=IPT1.NUM(3,IEL) 85 : DO I1=1,IDIM 86 : XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0* 87 : & (XCOOR((NN3-1)*(IDIM+1)+I1)+ 88 : & XCOOR((NN1-1)*(IDIM+1)+I1)) 89 : ENDDO 90 : ENDDO 91 : ELSEIF(IPT1.ITYPEL .EQ. 7)THEN 92 : C 93 : C********** TRI7 94 : C 95 : IF(LOTRI7)THEN 96 : C Elt already referred 97 : WRITE(IOIMP,*) 'Subroutine kdom2' 98 : WRITE(IOIMP,*) 'Mesh type not recognized' 98 : CALL ERREUR(5) 99 : ENDIF 100 : LOTRI7=.TRUE. 101 : C 102 : DO IEL=1,NBEL,1 103 : NN1=IPT1.NUM(1,IEL) 104 : NN2=IPT1.NUM(2,IEL) 105 : NN3=IPT1.NUM(3,IEL) 106 : NN4=IPT1.NUM(4,IEL) 107 : NN5=IPT1.NUM(5,IEL) 108 : NN6=IPT1.NUM(6,IEL) 109 : NN7=IPT1.NUM(7,IEL) 110 : DO I1=1,IDIM 111 : XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0* 112 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 113 : & XCOOR((NN3-1)*(IDIM+1)+I1)) 114 : XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0* 115 : & (XCOOR((NN3-1)*(IDIM+1)+I1)+ 116 : & XCOOR((NN5-1)*(IDIM+1)+I1)) 117 : XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0* 118 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 119 : & XCOOR((NN5-1)*(IDIM+1)+I1)) 120 : XCOOR((NN7-1)*(IDIM+1)+I1)= 121 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 122 : & XCOOR((NN3-1)*(IDIM+1)+I1)+ 123 : & XCOOR((NN5-1)*(IDIM+1)+I1))/3.0D0 124 : ENDDO 125 : ENDDO 126 : ELSEIF(IPT1.ITYPEL .EQ. 11)THEN 127 : C 128 : C********** QUA9 129 : C 130 : IF(LOQUA9)THEN 131 : C Elt already referred 132 : WRITE(IOIMP,*) 'Subroutine kdom2' 133 : WRITE(IOIMP,*) 'Mesh type not recognized' 133 : CALL ERREUR(5) 134 : ENDIF 135 : LOQUA9=.TRUE. 136 : C 137 : DO IEL=1,NBEL,1 138 : NN1=IPT1.NUM(1,IEL) 139 : NN2=IPT1.NUM(2,IEL) 140 : NN3=IPT1.NUM(3,IEL) 141 : NN4=IPT1.NUM(4,IEL) 142 : NN5=IPT1.NUM(5,IEL) 143 : NN6=IPT1.NUM(6,IEL) 144 : NN7=IPT1.NUM(7,IEL) 145 : NN8=IPT1.NUM(8,IEL) 146 : NN9=IPT1.NUM(9,IEL) 147 : C 148 : C**************** Centre de l'elt 149 : C 149 : CALL KDOM7(NN1,NN3,NN5,NN7,XCEN) 150 : IF(IERR.NE.0)GOTO 9999 151 : DO I1=1,IDIM 152 : XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0* 153 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 154 : & XCOOR((NN3-1)*(IDIM+1)+I1)) 155 : XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0* 156 : & (XCOOR((NN3-1)*(IDIM+1)+I1)+ 157 : & XCOOR((NN5-1)*(IDIM+1)+I1)) 158 : XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0* 159 : & (XCOOR((NN7-1)*(IDIM+1)+I1)+ 160 : & XCOOR((NN5-1)*(IDIM+1)+I1)) 161 : XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0* 162 : & (XCOOR((NN7-1)*(IDIM+1)+I1)+ 163 : & XCOOR((NN1-1)*(IDIM+1)+I1)) 164 : XCOOR((NN9-1)*(IDIM+1)+I1)=XCEN(I1) 165 : ENDDO 166 : ENDDO 167 : ELSEIF(IPT1.ITYPEL .EQ. 35)THEN 168 : C 169 : C********** TE15 170 : C 171 : IF(LOTE15)THEN 172 : C Elt already referred 173 : WRITE(IOIMP,*) 'Subroutine kdom2' 174 : WRITE(IOIMP,*) 'Mesh type not recognized' 174 : CALL ERREUR(5) 175 : ENDIF 176 : LOTE15=.TRUE. 177 : C 178 : DO IEL=1,NBEL,1 179 : NN1=IPT1.NUM(1,IEL) 180 : NN2=IPT1.NUM(2,IEL) 181 : NN3=IPT1.NUM(3,IEL) 182 : NN4=IPT1.NUM(4,IEL) 183 : NN5=IPT1.NUM(5,IEL) 184 : NN6=IPT1.NUM(6,IEL) 185 : NN7=IPT1.NUM(7,IEL) 186 : NN8=IPT1.NUM(8,IEL) 187 : NN9=IPT1.NUM(9,IEL) 188 : NN10=IPT1.NUM(10,IEL) 189 : NN11=IPT1.NUM(11,IEL) 190 : NN12=IPT1.NUM(12,IEL) 191 : NN13=IPT1.NUM(13,IEL) 192 : NN14=IPT1.NUM(14,IEL) 193 : NN15=IPT1.NUM(15,IEL) 194 : DO I1=1,IDIM 195 : XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0* 196 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 197 : & XCOOR((NN3-1)*(IDIM+1)+I1)) 198 : XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0* 199 : & (XCOOR((NN3-1)*(IDIM+1)+I1)+ 200 : & XCOOR((NN5-1)*(IDIM+1)+I1)) 201 : XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0* 202 : & (XCOOR((NN5-1)*(IDIM+1)+I1)+ 203 : & XCOOR((NN1-1)*(IDIM+1)+I1)) 204 : XCOOR((NN7-1)*(IDIM+1)+I1)=0.5D0* 205 : & (XCOOR((NN10-1)*(IDIM+1)+I1)+ 206 : & XCOOR((NN1-1)*(IDIM+1)+I1)) 207 : XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0* 208 : & (XCOOR((NN10-1)*(IDIM+1)+I1)+ 209 : & XCOOR((NN3-1)*(IDIM+1)+I1)) 210 : XCOOR((NN9-1)*(IDIM+1)+I1)=0.5D0* 211 : & (XCOOR((NN10-1)*(IDIM+1)+I1)+ 212 : & XCOOR((NN5-1)*(IDIM+1)+I1)) 213 : XCOOR((NN11-1)*(IDIM+1)+I1)= 214 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 215 : & XCOOR((NN3-1)*(IDIM+1)+I1)+ 216 : & XCOOR((NN5-1)*(IDIM+1)+I1))/3.0D0 217 : XCOOR((NN12-1)*(IDIM+1)+I1)= 218 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 219 : & XCOOR((NN3-1)*(IDIM+1)+I1)+ 220 : & XCOOR((NN10-1)*(IDIM+1)+I1))/3.0D0 221 : XCOOR((NN13-1)*(IDIM+1)+I1)= 222 : & (XCOOR((NN3-1)*(IDIM+1)+I1)+ 223 : & XCOOR((NN5-1)*(IDIM+1)+I1)+ 224 : & XCOOR((NN10-1)*(IDIM+1)+I1))/3.0D0 225 : XCOOR((NN14-1)*(IDIM+1)+I1)= 226 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 227 : & XCOOR((NN5-1)*(IDIM+1)+I1)+ 228 : & XCOOR((NN10-1)*(IDIM+1)+I1))/3.0D0 229 : XCOOR((NN15-1)*(IDIM+1)+I1)= 230 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 231 : & XCOOR((NN3-1)*(IDIM+1)+I1)+ 232 : & XCOOR((NN5-1)*(IDIM+1)+I1)+ 233 : & XCOOR((NN10-1)*(IDIM+1)+I1))/4.0D0 234 : ENDDO 235 : ENDDO 236 : ELSEIF(IPT1.ITYPEL .EQ. 36)THEN 237 : C 238 : C********** PY19 239 : C 240 : IF(LOPY19)THEN 241 : C Elt already referred 242 : WRITE(IOIMP,*) 'Subroutine kdom2' 243 : WRITE(IOIMP,*) 'Mesh type not recognized' 243 : CALL ERREUR(5) 244 : ENDIF 245 : LOPY19=.TRUE. 246 : C 247 : DO IEL=1,NBEL,1 248 : NN1=IPT1.NUM(1,IEL) 249 : NN2=IPT1.NUM(2,IEL) 250 : NN3=IPT1.NUM(3,IEL) 251 : NN4=IPT1.NUM(4,IEL) 252 : NN5=IPT1.NUM(5,IEL) 253 : NN6=IPT1.NUM(6,IEL) 254 : NN7=IPT1.NUM(7,IEL) 255 : NN8=IPT1.NUM(8,IEL) 256 : NN9=IPT1.NUM(9,IEL) 257 : NN10=IPT1.NUM(10,IEL) 258 : NN11=IPT1.NUM(11,IEL) 259 : NN12=IPT1.NUM(12,IEL) 260 : NN13=IPT1.NUM(13,IEL) 261 : NN14=IPT1.NUM(14,IEL) 262 : NN15=IPT1.NUM(15,IEL) 263 : NN16=IPT1.NUM(16,IEL) 264 : NN17=IPT1.NUM(17,IEL) 265 : NN18=IPT1.NUM(18,IEL) 266 : NN19=IPT1.NUM(19,IEL) 267 : C Computation of the center position 267 : CALL KDOM7(NN1,NN3,NN5,NN7,XCEN) 268 : IF(IERR.NE.0)GOTO 9999 269 : DO I1=1,IDIM 270 : XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0* 271 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 272 : & XCOOR((NN3-1)*(IDIM+1)+I1)) 273 : XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0* 274 : & (XCOOR((NN3-1)*(IDIM+1)+I1)+ 275 : & XCOOR((NN5-1)*(IDIM+1)+I1)) 276 : XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0* 277 : & (XCOOR((NN5-1)*(IDIM+1)+I1)+ 278 : & XCOOR((NN7-1)*(IDIM+1)+I1)) 279 : XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0* 280 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 281 : & XCOOR((NN7-1)*(IDIM+1)+I1)) 282 : XCOOR((NN9-1)*(IDIM+1)+I1)=0.5D0* 283 : & (XCOOR((NN13-1)*(IDIM+1)+I1)+ 284 : & XCOOR((NN1-1)*(IDIM+1)+I1)) 285 : XCOOR((NN10-1)*(IDIM+1)+I1)=0.5D0* 286 : & (XCOOR((NN13-1)*(IDIM+1)+I1)+ 287 : & XCOOR((NN3-1)*(IDIM+1)+I1)) 288 : XCOOR((NN11-1)*(IDIM+1)+I1)=0.5D0* 289 : & (XCOOR((NN13-1)*(IDIM+1)+I1)+ 290 : & XCOOR((NN5-1)*(IDIM+1)+I1)) 291 : XCOOR((NN12-1)*(IDIM+1)+I1)=0.5D0* 292 : & (XCOOR((NN13-1)*(IDIM+1)+I1)+ 293 : & XCOOR((NN7-1)*(IDIM+1)+I1)) 294 : XCOOR((NN14-1)*(IDIM+1)+I1)=XCEN(I1) 295 : XCOOR((NN15-1)*(IDIM+1)+I1)= 296 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 297 : & XCOOR((NN3-1)*(IDIM+1)+I1)+ 298 : & XCOOR((NN13-1)*(IDIM+1)+I1))/3.0D0 299 : XCOOR((NN16-1)*(IDIM+1)+I1)= 300 : & (XCOOR((NN3-1)*(IDIM+1)+I1)+ 301 : & XCOOR((NN5-1)*(IDIM+1)+I1)+ 302 : & XCOOR((NN13-1)*(IDIM+1)+I1))/3.0D0 303 : XCOOR((NN17-1)*(IDIM+1)+I1)= 304 : & (XCOOR((NN5-1)*(IDIM+1)+I1)+ 305 : & XCOOR((NN7-1)*(IDIM+1)+I1)+ 306 : & XCOOR((NN13-1)*(IDIM+1)+I1))/3.0D0 307 : XCOOR((NN18-1)*(IDIM+1)+I1)= 308 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 309 : & XCOOR((NN7-1)*(IDIM+1)+I1)+ 310 : & XCOOR((NN13-1)*(IDIM+1)+I1))/3.0D0 311 : ENDDO 312 : C 313 : C************* To compute the center, we divide the pyramid into 314 : C 4 tetrahedra (see kdom3) 315 : C 315 : CALL KDOM3(NN1,NN3,NN5,NN7,NN14,NN13,XCEN,VOL) 316 : DO I1=1,3,1 317 : XCOOR((NN19-1)*(IDIM+1)+I1)=XCEN(I1) 318 : ENDDO 319 : ENDDO 320 : ELSEIF(IPT1.ITYPEL .EQ. 34)THEN 321 : C 322 : C********** PR21 323 : C 324 : IF(LOPR21)THEN 325 : C Elt already referred 326 : WRITE(IOIMP,*) 'Subroutine kdom2' 327 : WRITE(IOIMP,*) 'Mesh type not recognized' 327 : CALL ERREUR(5) 328 : ENDIF 329 : LOPR21=.TRUE. 330 : C 331 : DO IEL=1,NBEL,1 332 : NN1=IPT1.NUM(1,IEL) 333 : NN2=IPT1.NUM(2,IEL) 334 : NN3=IPT1.NUM(3,IEL) 335 : NN4=IPT1.NUM(4,IEL) 336 : NN5=IPT1.NUM(5,IEL) 337 : NN6=IPT1.NUM(6,IEL) 338 : NN7=IPT1.NUM(7,IEL) 339 : NN8=IPT1.NUM(8,IEL) 340 : NN9=IPT1.NUM(9,IEL) 341 : NN10=IPT1.NUM(10,IEL) 342 : NN11=IPT1.NUM(11,IEL) 343 : NN12=IPT1.NUM(12,IEL) 344 : NN13=IPT1.NUM(13,IEL) 345 : NN14=IPT1.NUM(14,IEL) 346 : NN15=IPT1.NUM(15,IEL) 347 : NN16=IPT1.NUM(16,IEL) 348 : NN17=IPT1.NUM(17,IEL) 349 : NN18=IPT1.NUM(18,IEL) 350 : NN19=IPT1.NUM(19,IEL) 351 : NN20=IPT1.NUM(20,IEL) 352 : NN21=IPT1.NUM(21,IEL) 353 : DO I1=1,IDIM 354 : XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0* 355 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 356 : & XCOOR((NN3-1)*(IDIM+1)+I1)) 357 : XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0* 358 : & (XCOOR((NN3-1)*(IDIM+1)+I1)+ 359 : & XCOOR((NN5-1)*(IDIM+1)+I1)) 360 : XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0* 361 : & (XCOOR((NN5-1)*(IDIM+1)+I1)+ 362 : & XCOOR((NN1-1)*(IDIM+1)+I1)) 363 : XCOOR((NN7-1)*(IDIM+1)+I1)=0.5D0* 364 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 365 : & XCOOR((NN10-1)*(IDIM+1)+I1)) 366 : XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0* 367 : & (XCOOR((NN3-1)*(IDIM+1)+I1)+ 368 : & XCOOR((NN12-1)*(IDIM+1)+I1)) 369 : XCOOR((NN9-1)*(IDIM+1)+I1)=0.5D0* 370 : & (XCOOR((NN5-1)*(IDIM+1)+I1)+ 371 : & XCOOR((NN14-1)*(IDIM+1)+I1)) 372 : XCOOR((NN11-1)*(IDIM+1)+I1)=0.5D0* 373 : & (XCOOR((NN10-1)*(IDIM+1)+I1)+ 374 : & XCOOR((NN12-1)*(IDIM+1)+I1)) 375 : XCOOR((NN13-1)*(IDIM+1)+I1)=0.5D0* 376 : & (XCOOR((NN12-1)*(IDIM+1)+I1)+ 377 : & XCOOR((NN14-1)*(IDIM+1)+I1)) 378 : XCOOR((NN15-1)*(IDIM+1)+I1)=0.5D0* 379 : & (XCOOR((NN10-1)*(IDIM+1)+I1)+ 380 : & XCOOR((NN14-1)*(IDIM+1)+I1)) 381 : XCOOR((NN19-1)*(IDIM+1)+I1)= 382 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 383 : & XCOOR((NN3-1)*(IDIM+1)+I1)+ 384 : & XCOOR((NN5-1)*(IDIM+1)+I1))/3.0D0 385 : XCOOR((NN20-1)*(IDIM+1)+I1)= 386 : & (XCOOR((NN10-1)*(IDIM+1)+I1)+ 387 : & XCOOR((NN12-1)*(IDIM+1)+I1)+ 388 : & XCOOR((NN14-1)*(IDIM+1)+I1))/3.0D0 389 : ENDDO 390 : C Computation of the QUA center positions 390 : CALL KDOM7(NN1,NN10,NN12,NN3,XCEN) 391 : IF(IERR.NE.0)GOTO 9999 392 : DO I1=1,IDIM,1 393 : XCOOR((NN16-1)*(IDIM+1)+I1)=XCEN(I1) 394 : ENDDO 394 : CALL KDOM7(NN5,NN3,NN12,NN14,XCEN) 395 : IF(IERR.NE.0)GOTO 9999 396 : DO I1=1,IDIM,1 397 : XCOOR((NN17-1)*(IDIM+1)+I1)=XCEN(I1) 398 : ENDDO 398 : CALL KDOM7(NN10,NN1,NN5,NN14,XCEN) 399 : IF(IERR.NE.0)GOTO 9999 400 : DO I1=1,IDIM,1 401 : XCOOR((NN18-1)*(IDIM+1)+I1)=XCEN(I1) 402 : ENDDO 403 : C 404 : C************* To compute the center, we divide the PR21 405 : C into 2 tetrahedra and 3 pyramids 406 : C (1,3,5,PC), (10,14,12,PC), (1,10,12,3,16,PC), 407 : C (3,12,14,5,17,PC),(1,5,14,10,18,PC) 408 : C | 409 : C center of 3,12,14,5 410 : C The results does not depend of the choice of 411 : C PC. 412 : C 413 : C 414 : C 415 : DO I1=1,3,1 416 : XCOOR((NN21-1)*(IDIM+1)+I1)= 417 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 418 : & XCOOR((NN3-1)*(IDIM+1)+I1)+ 419 : & XCOOR((NN5-1)*(IDIM+1)+I1)+ 420 : & XCOOR((NN10-1)*(IDIM+1)+I1)+ 421 : & XCOOR((NN12-1)*(IDIM+1)+I1)+ 422 : & XCOOR((NN14-1)*(IDIM+1)+I1))/6.0D0 423 : 424 : ENDDO 425 : C 425 : CALL KDOM3(NN10,NN12,NN3,NN1,NN16,NN21,X1,VOL1) 425 : CALL KDOM3(NN3,NN12,NN14,NN5,NN17,NN21,X2,VOL2) 425 : CALL KDOM3(NN1,NN5,NN14,NN10,NN18,NN21,X3,VOL3) 425 : CALL KDOM4(NN1,NN3,NN5,NN21,X4,VOL4) 425 : CALL KDOM4(NN10,NN14,NN12,NN21,X5,VOL5) 426 : VOL=VOL1+VOL2+VOL3+VOL4+VOL5 427 : C 428 : DO I1=1,3,1 429 : XCOOR((NN21-1)*(IDIM+1)+I1)=(VOL1*X1(I1)+VOL2*X2(I1) 430 : & +VOL3*X3(I1)+VOL4*X4(I1)+VOL5*X5(I1))/VOL 431 : ENDDO 432 : C 433 : C********** We check whether the center of gravity of the PR21 is 434 : C inside its domain 435 : C 435 : CALL KDOM8(NN1,NN2,NN3,NN4,NN5,NN6,NN7,NN8,NN9,NN10, 436 : & NN11,NN12,NN13,NN14,NN15,NN16,NN17,NN18,NN19,NN20, 437 : & NN21) 438 : IF(IERR.NE.0)GOTO 9999 439 : ENDDO 440 : IF(NBSOUS .NE. 1) SEGDES IPT1 441 : ELSEIF(IPT1.ITYPEL .EQ. 33)THEN 442 : C 443 : C********** CU27 444 : C 445 : IF(LOCU27)THEN 446 : C Elt already referred 447 : WRITE(IOIMP,*) 'Subroutine kdom2' 448 : WRITE(IOIMP,*) 'Mesh type not recognized' 448 : CALL ERREUR(5) 449 : ENDIF 450 : LOCU27=.TRUE. 451 : C 452 : DO IEL=1,NBEL,1 453 : NN1=IPT1.NUM(1,IEL) 454 : NN2=IPT1.NUM(2,IEL) 455 : NN3=IPT1.NUM(3,IEL) 456 : NN4=IPT1.NUM(4,IEL) 457 : NN5=IPT1.NUM(5,IEL) 458 : NN6=IPT1.NUM(6,IEL) 459 : NN7=IPT1.NUM(7,IEL) 460 : NN8=IPT1.NUM(8,IEL) 461 : NN9=IPT1.NUM(9,IEL) 462 : NN10=IPT1.NUM(10,IEL) 463 : NN11=IPT1.NUM(11,IEL) 464 : NN12=IPT1.NUM(12,IEL) 465 : NN13=IPT1.NUM(13,IEL) 466 : NN14=IPT1.NUM(14,IEL) 467 : NN15=IPT1.NUM(15,IEL) 468 : NN16=IPT1.NUM(16,IEL) 469 : NN17=IPT1.NUM(17,IEL) 470 : NN18=IPT1.NUM(18,IEL) 471 : NN19=IPT1.NUM(19,IEL) 472 : NN20=IPT1.NUM(20,IEL) 473 : NN21=IPT1.NUM(21,IEL) 474 : NN22=IPT1.NUM(22,IEL) 475 : NN23=IPT1.NUM(23,IEL) 476 : NN24=IPT1.NUM(24,IEL) 477 : NN25=IPT1.NUM(25,IEL) 478 : NN26=IPT1.NUM(26,IEL) 479 : NN27=IPT1.NUM(27,IEL) 480 : DO I1=1,IDIM 481 : XCOOR((NN2-1)*(IDIM+1)+I1)=0.5D0* 482 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 483 : & XCOOR((NN3-1)*(IDIM+1)+I1)) 484 : XCOOR((NN4-1)*(IDIM+1)+I1)=0.5D0* 485 : & (XCOOR((NN3-1)*(IDIM+1)+I1)+ 486 : & XCOOR((NN5-1)*(IDIM+1)+I1)) 487 : XCOOR((NN6-1)*(IDIM+1)+I1)=0.5D0* 488 : & (XCOOR((NN5-1)*(IDIM+1)+I1)+ 489 : & XCOOR((NN7-1)*(IDIM+1)+I1)) 490 : XCOOR((NN8-1)*(IDIM+1)+I1)=0.5D0* 491 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 492 : & XCOOR((NN7-1)*(IDIM+1)+I1)) 493 : XCOOR((NN9-1)*(IDIM+1)+I1)=0.5D0* 494 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 495 : & XCOOR((NN13-1)*(IDIM+1)+I1)) 496 : XCOOR((NN10-1)*(IDIM+1)+I1)=0.5D0* 497 : & (XCOOR((NN3-1)*(IDIM+1)+I1)+ 498 : & XCOOR((NN15-1)*(IDIM+1)+I1)) 499 : XCOOR((NN11-1)*(IDIM+1)+I1)=0.5D0* 500 : & (XCOOR((NN5-1)*(IDIM+1)+I1)+ 501 : & XCOOR((NN17-1)*(IDIM+1)+I1)) 502 : XCOOR((NN12-1)*(IDIM+1)+I1)=0.5D0* 503 : & (XCOOR((NN7-1)*(IDIM+1)+I1)+ 504 : & XCOOR((NN19-1)*(IDIM+1)+I1)) 505 : XCOOR((NN14-1)*(IDIM+1)+I1)=0.5D0* 506 : & (XCOOR((NN13-1)*(IDIM+1)+I1)+ 507 : & XCOOR((NN15-1)*(IDIM+1)+I1)) 508 : XCOOR((NN16-1)*(IDIM+1)+I1)=0.5D0* 509 : & (XCOOR((NN17-1)*(IDIM+1)+I1)+ 510 : & XCOOR((NN15-1)*(IDIM+1)+I1)) 511 : XCOOR((NN18-1)*(IDIM+1)+I1)=0.5D0* 512 : & (XCOOR((NN17-1)*(IDIM+1)+I1)+ 513 : & XCOOR((NN19-1)*(IDIM+1)+I1)) 514 : XCOOR((NN20-1)*(IDIM+1)+I1)=0.5D0* 515 : & (XCOOR((NN19-1)*(IDIM+1)+I1)+ 516 : & XCOOR((NN13-1)*(IDIM+1)+I1)) 517 : ENDDO 518 : C Computation of the QUA center positions 518 : CALL KDOM7(NN1,NN13,NN15,NN3,XCEN) 519 : IF(IERR.NE.0) GOTO 9999 520 : DO I1=1,IDIM,1 521 : XCOOR((NN21-1)*(IDIM+1)+I1)=XCEN(I1) 522 : ENDDO 522 : CALL KDOM7(NN3,NN15,NN17,NN5,XCEN) 523 : IF(IERR.NE.0)GOTO 9999 524 : DO I1=1,IDIM,1 525 : XCOOR((NN22-1)*(IDIM+1)+I1)=XCEN(I1) 526 : ENDDO 526 : CALL KDOM7(NN5,NN17,NN19,NN7,XCEN) 527 : IF(IERR.NE.0)GOTO 9999 528 : DO I1=1,IDIM,1 529 : XCOOR((NN23-1)*(IDIM+1)+I1)=XCEN(I1) 530 : ENDDO 530 : CALL KDOM7(NN7,NN19,NN13,NN1,XCEN) 531 : IF(IERR.NE.0)GOTO 9999 532 : DO I1=1,IDIM,1 533 : XCOOR((NN24-1)*(IDIM+1)+I1)=XCEN(I1) 534 : ENDDO 534 : CALL KDOM7(NN1,NN3,NN5,NN7,XCEN) 535 : IF(IERR.NE.0)GOTO 9999 536 : DO I1=1,IDIM,1 537 : XCOOR((NN25-1)*(IDIM+1)+I1)=XCEN(I1) 538 : ENDDO 538 : CALL KDOM7(NN13,NN19,NN17,NN15,XCEN) 539 : IF(IERR.NE.0)GOTO 9999 540 : DO I1=1,IDIM,1 541 : XCOOR((NN26-1)*(IDIM+1)+I1)=XCEN(I1) 542 : ENDDO 543 : C 544 : C************* To compute the center, we divide the cube 545 : C into 6 pyramids. 546 : C (1,3,5,7,25,PC) 547 : C (13,19,17,15,26,PC) 548 : C (1,13,15,3,21,PC) 549 : C (3,15,17,5,22,PC) 550 : C (5,17,19,7,23,PC) 551 : C (1,7,19,13,24,PC) 552 : 553 : C The results does not depend of the choice of 554 : C PC. 555 : C 556 : DO I1=1,3,1 557 : XCOOR((NN27-1)*(IDIM+1)+I1)= 558 : & (XCOOR((NN1-1)*(IDIM+1)+I1)+ 559 : & XCOOR((NN3-1)*(IDIM+1)+I1)+ 560 : & XCOOR((NN5-1)*(IDIM+1)+I1)+ 561 : & XCOOR((NN7-1)*(IDIM+1)+I1)+ 562 : & XCOOR((NN13-1)*(IDIM+1)+I1)+ 563 : & XCOOR((NN15-1)*(IDIM+1)+I1)+ 564 : & XCOOR((NN17-1)*(IDIM+1)+I1)+ 565 : & XCOOR((NN19-1)*(IDIM+1)+I1))/8.0D0 566 : 567 : ENDDO 568 : C 568 : CALL KDOM3(NN1,NN3,NN5,NN7,NN25,NN27,X1,VOL1) 568 : CALL KDOM3(NN13,NN19,NN17,NN15,NN26,NN27,X2,VOL2) 568 : CALL KDOM3(NN1,NN13,NN15,NN3,NN21,NN27,X3,VOL3) 568 : CALL KDOM3(NN3,NN15,NN17,NN5,NN22,NN27,X4,VOL4) 568 : CALL KDOM3(NN5,NN17,NN19,NN7,NN23,NN27,X5,VOL5) 568 : CALL KDOM3(NN1,NN7,NN19,NN13,NN24,NN27,X6,VOL6) 569 : VOL=VOL1+VOL2+VOL3+VOL4+VOL5+VOL6 570 : DO I1=1,3,1 571 : XCOOR((NN27-1)*(IDIM+1)+I1)=(VOL1*X1(I1)+VOL2*X2(I1) 572 : & +VOL3*X3(I1)+VOL4*X4(I1)+VOL5*X5(I1)+VOL6*X6(I1)) 573 : & /VOL 574 : ENDDO 575 : C 576 : C********** We check whether the center of gravity of the CU27 is 577 : C inside its domain 578 : C 578 : CALL KDOM9(NN1,NN2,NN3,NN4,NN5,NN6,NN7,NN8,NN9,NN10, 579 : & NN11,NN12,NN13,NN14,NN15,NN16,NN17,NN18,NN19,NN20, 580 : & NN21,NN22,NN23,NN24,NN25,NN26,NN27) 581 : IF(IERR.NE.0)GOTO 9999 582 : 583 : ENDDO 584 : ELSE 585 : WRITE(IOIMP,*) 'Euler model' 585 : CALL ERREUR(21) 586 : ENDIF 587 : IF(NBSOUS .NE. 1) SEGDES IPT1 588 : ENDDO 589 : C 590 : SEGDES MELEME 591 : C 592 : 9999 RETURN 593 : C 594 : END 595 : 596 :
© Cast3M 2003 - All rights reserved.
Disclaimer