1 : C FACED2 SOURCE PV 12/04/13 21:15:20 7351 2 : C TRACE D'ISOVALEUR EN COMMENCANT PAR CELLES DE DERRIERE 3 : C 4 : SUBROUTINE FACED2(MELEME,XPROJ,ICPR,VCHC,VCPCHA,PTI,NISO,IVU, 5 : # MCOUP,mcham) 6 : * si mchaml est non nul on trace un chamelem aux noeuds 7 : * sinon c'est un champoint 8 : IMPLICIT INTEGER(I-N) 9 : -INC CCOPTIO 10 : -INC SMELEME 11 : -INC CCGEOME 12 : -INC SMCHAML 13 : DIMENSION XX(13),YY(13),ZZ(13),VV(13) 14 : SEGMENT XPROJ(3,1) 15 : SEGMENT ICPR(1) 16 : SEGMENT VCPCHA(XCOOR(/1)/(IDIM+1)) 17 : SEGMENT ISEG2(2,NSEG2) 18 : SEGMENT ISEG3(3,NSEG3) 19 : SEGMENT IFAC3(4,NFAC3) 20 : SEGMENT IFAC4(5,NFAC4) 21 : SEGMENT IFAC6(7,NFAC6) 22 : SEGMENT IFAC8(9,NFAC8) 23 : SEGMENT IFACP(15,NFACP) 24 : SEGMENT XSEG2(2,NSEG2) 25 : SEGMENT XSEG3(3,NSEG3) 26 : SEGMENT XFAC3(3,NFAC3) 27 : SEGMENT XFAC4(4,NFAC4) 28 : SEGMENT XFAC6(6,NFAC6) 29 : SEGMENT XFAC8(8,NFAC8) 30 : SEGMENT XFACP(14,NFACP) 31 : SEGMENT NSOMP(NFACP) 32 : SEGMENT TFAC(NFAC) 33 : SEGMENT KFAK(NFAC) 34 : SEGMENT NAUX(max(2,NFAC)) 35 : SEGMENT IVU(0) 36 : SEGMENT MCOUP(0) 37 : DIMENSION XTR(8),YTR(8),ZTR(8) 38 : REAL VCHC 39 : DIMENSION VCHC(*) 40 : 41 : SEGACT MELEME 42 : * ipt1=lisous(1) 43 : * segact ipt1 44 : * SEGACT XPROJ,ICPR 45 : MELSAU=MELEME 46 : IPT1=MELEME 47 : melval=ipt1 48 : * write (6,*) ' faced2 velche ' 49 : * write (6,*) (velche(ix,1),ix=1,20) 50 : if (mcham.eq.0) then 50 : CALL ECROBJ('MAILLAGE',IPT1) 50 : CALL ENVELO 50 : CALL LIROBJ('MAILLAGE',MELEME,1,IRETOU) 51 : else 51 : CALL ENVEL1(ipt1,meleme) 52 : endif 53 : IF (IERR.NE.0) RETURN 54 : 55 : SEGACT MELEME 56 : SEGACT XPROJ,ICPR 57 : NBPOIN=XPROJ(/2) 58 : TMIN=1E30 59 : TMAX=-1E30 60 : DO 1 I=1,NBPOIN 61 : TMIN=MIN(TMIN,XPROJ(3,I)) 62 : TMAX=MAX(TMAX,XPROJ(3,I)) 63 : 1 CONTINUE 64 : TDIST=TMAX-TMIN 65 : NSEG2=0 66 : NSEG3=0 67 : NFAC3=0 68 : NFAC4=0 69 : NFAC6=0 70 : NFAC8=0 71 : NFACP=0 72 : IPT1=MELEME 73 : SEGACT MELEME 74 : DO 10 IOB=1,MAX(1,LISOUS(/1)) 75 : IF (LISOUS(/1).NE.0) THEN 76 : IPT1=LISOUS(IOB) 77 : SEGACT IPT1 78 : ENDIF 79 : NBELEM=IPT1.NUM(/2) 80 : ILTEL=LTEL(1,IPT1.ITYPEL) 81 : IF (ILTEL.EQ.0) GOTO 12 82 : ILTAD=LTEL(2,IPT1.ITYPEL) 83 : DO 13 IF=1,ILTEL 84 : IFT=LDEL(1,ILTAD+IF-1) 85 : GOTO (21,22,23,24,25,26),IFT 86 : 21 NFAC3=NFAC3+NBELEM 87 : GOTO 27 88 : 22 NFAC4=NFAC4+NBELEM 89 : GOTO 27 90 : 23 NFAC6=NFAC6+NBELEM 91 : GOTO 27 92 : 24 NFAC8=NFAC8+NBELEM 93 : GOTO 27 94 : 25 GOTO 27 95 : 26 NFACP=NFACP+NBELEM 96 : GOTO 27 97 : 27 CONTINUE 98 : 13 CONTINUE 99 : goto 10 100 : 12 CONTINUE 101 : if (ipt1.itypel.eq.2) nseg2=nseg2+nbelem 102 : if (ipt1.itypel.eq.3) nseg3=nseg3+nbelem 103 : 10 CONTINUE 104 : SEGINI iseg2,iseg3,IFAC3,IFAC4,IFAC6,IFAC8,IFACP,NSOMP 105 : if (mcham.ne.0) SEGINI xseg2,xseg3,xFAC3,xFAC4,xFAC6,xFAC8,xFACP 106 : NSEG2=0 107 : NSEG3=0 108 : NFAC3=0 109 : NFAC4=0 110 : NFAC6=0 111 : NFAC8=0 112 : NFACP=0 113 : ICOUPE=0 114 : DO 50 IOB=1,MAX(1,LISOUS(/1)) 115 : IF (LISOUS(/1).NE.0) THEN 116 : IPT1=LISOUS(IOB) 117 : ivm1=0 118 : ivm2=0 119 : melval=meleme 120 : if (mcham.ne.0) then 121 : melval=lisref(iob) 122 : ivm1=velche(/1) 123 : ivm2=velche(/2) 124 : if (ivm1*ivm2.eq.0) return 125 : endif 126 : ICOUPE=0 127 : IF (IOB.EQ.LISOUS(/1).AND.MCOUP.NE.0) ICOUPE=1 128 : ENDIF 129 : NBELEM=IPT1.NUM(/2) 130 : NBNN =IPT1.NUM(/1) 131 : ILTEL=LTEL(1,IPT1.ITYPEL) 132 : IF (ILTEL.EQ.0) GOTO 52 133 : ILTAD=LTEL(2,IPT1.ITYPEL) 134 : DO 60 IF=1,ILTEL 135 : ITYP=LDEL(1,ILTAD+IF-1) 136 : IAD=LDEL(2,ILTAD+IF-1) 137 : DO 80 IEL=1,NBELEM 138 : do 70 inn=1,nbnn 139 : if (icpr(ipt1.num(inn,iel)).eq.0) then 139 : call erreur(23) 140 : goto 80 141 : endif 142 : 70 continue 143 : GOTO (61,62,63,64,67,66),ITYP 144 : 61 NFAC3=NFAC3+1 145 : IFAC3(1,NFAC3)=IPT1.NUM(LFAC(IAD),IEL) 146 : IFAC3(2,NFAC3)=IPT1.NUM(LFAC(IAD+1),IEL) 147 : IFAC3(3,NFAC3)=IPT1.NUM(LFAC(IAD+2),IEL) 148 : if (mcham.ne.0) then 149 : xFAC3(1,NFAC3)=velche(min(ivm1,LFAC(IAD)),min(ivm2,IEL)) 150 : xFAC3(2,NFAC3)=velche(min(ivm1,LFAC(IAD+1)),min(ivm2,IEL)) 151 : xFAC3(3,NFAC3)=velche(min(ivm1,LFAC(IAD+2)),min(ivm2,IEL)) 152 : endif 153 : IFAC3(4,NFAC3)=0 154 : IF (IVU(ICPR(IFAC3(1,NFAC3))).EQ.1.AND. 155 : # IVU(ICPR(IFAC3(2,NFAC3))).EQ.1.AND. 156 : # IVU(ICPR(IFAC3(3,NFAC3))).EQ.1) IFAC3(4,NFAC3)=1 157 : IF (ICOUPE.EQ.1) THEN 158 : IF (MCOUP(IEL)/8.EQ.1) IFAC3(4,NFAC3)=2 159 : IF (MCOUP(IEL)/16.EQ.1) IFAC3(4,NFAC3)=3 160 : ENDIF 161 : GOTO 67 162 : 62 NFAC4=NFAC4+1 163 : IFAC4(1,NFAC4)=IPT1.NUM(LFAC(IAD),IEL) 164 : IFAC4(2,NFAC4)=IPT1.NUM(LFAC(IAD+1),IEL) 165 : IFAC4(3,NFAC4)=IPT1.NUM(LFAC(IAD+2),IEL) 166 : IFAC4(4,NFAC4)=IPT1.NUM(LFAC(IAD+3),IEL) 167 : if (mcham.ne.0) then 168 : xFAC4(1,NFAC4)=velche(min(ivm1,LFAC(IAD)),min(ivm2,IEL)) 169 : xFAC4(2,NFAC4)=velche(min(ivm1,LFAC(IAD+1)),min(ivm2,IEL)) 170 : xFAC4(3,NFAC4)=velche(min(ivm1,LFAC(IAD+2)),min(ivm2,IEL)) 171 : xFAC4(4,NFAC4)=velche(min(ivm1,LFAC(IAD+3)),min(ivm2,IEL)) 172 : endif 173 : IFAC4(5,NFAC4)=0 174 : IF (IVU(ICPR(IFAC4(1,NFAC4))).EQ.1.AND. 175 : # IVU(ICPR(IFAC4(2,NFAC4))).EQ.1.AND. 176 : # IVU(ICPR(IFAC4(3,NFAC4))).EQ.1.AND. 177 : # IVU(ICPR(IFAC4(4,NFAC4))).EQ.1) 178 : # IFAC4(5,NFAC4)=1 179 : GOTO 67 180 : 63 NFAC6=NFAC6+1 181 : IFAC6(1,NFAC6)=IPT1.NUM(LFAC(IAD) ,IEL) 182 : IFAC6(2,NFAC6)=IPT1.NUM(LFAC(IAD+1),IEL) 183 : IFAC6(3,NFAC6)=IPT1.NUM(LFAC(IAD+2),IEL) 184 : IFAC6(4,NFAC6)=IPT1.NUM(LFAC(IAD+3),IEL) 185 : IFAC6(5,NFAC6)=IPT1.NUM(LFAC(IAD+4),IEL) 186 : IFAC6(6,NFAC6)=IPT1.NUM(LFAC(IAD+5),IEL) 187 : if (mcham.ne.0) then 188 : xFAC6(1,NFAC6)=velche(min(ivm1,LFAC(IAD)) ,min(ivm2,IEL)) 189 : xFAC6(2,NFAC6)=velche(min(ivm1,LFAC(IAD+1)),min(ivm2,IEL)) 190 : xFAC6(3,NFAC6)=velche(min(ivm1,LFAC(IAD+2)),min(ivm2,IEL)) 191 : xFAC6(4,NFAC6)=velche(min(ivm1,LFAC(IAD+3)),min(ivm2,IEL)) 192 : xFAC6(5,NFAC6)=velche(min(ivm1,LFAC(IAD+4)),min(ivm2,IEL)) 193 : xFAC6(6,NFAC6)=velche(min(ivm1,LFAC(IAD+5)),min(ivm2,IEL)) 194 : endif 195 : IFAC6(7,NFAC6)=0 196 : ipa=icpr(ifac6(1,nfac6)) 197 : ipb=icpr(ifac6(2,nfac6)) 198 : ipc=icpr(ifac6(3,nfac6)) 199 : ipd=icpr(ifac6(4,nfac6)) 200 : ipe=icpr(ifac6(5,nfac6)) 201 : ipf=icpr(ifac6(6,nfac6)) 202 : IF (IVU(ipa).EQ.1.AND.IVU(ipb).EQ.1.AND.IVU(ipc).EQ.1 203 : > .AND.IVU(ipd).EQ.1.AND.IVU(ipe).EQ.1 204 : > .AND.IVU(ipf).EQ.1) 205 : > IFAC6(7,NFAC6)=1 206 : GOTO 67 207 : 64 CONTINUE 208 : NFAC8=NFAC8+1 209 : IFAC8(1,NFAC8)=IPT1.NUM(LFAC(IAD) ,IEL) 210 : IFAC8(2,NFAC8)=IPT1.NUM(LFAC(IAD+1),IEL) 211 : IFAC8(3,NFAC8)=IPT1.NUM(LFAC(IAD+2),IEL) 212 : IFAC8(4,NFAC8)=IPT1.NUM(LFAC(IAD+3),IEL) 213 : IFAC8(5,NFAC8)=IPT1.NUM(LFAC(IAD+4),IEL) 214 : IFAC8(6,NFAC8)=IPT1.NUM(LFAC(IAD+5),IEL) 215 : IFAC8(7,NFAC8)=IPT1.NUM(LFAC(IAD+6),IEL) 216 : IFAC8(8,NFAC8)=IPT1.NUM(LFAC(IAD+7),IEL) 217 : if (mcham.ne.0) then 218 : xFAC8(1,NFAC8)=velche(min(ivm1,LFAC(IAD)) ,min(ivm2,IEL)) 219 : xFAC8(2,NFAC8)=velche(min(ivm1,LFAC(IAD+1)),min(ivm2,IEL)) 220 : xFAC8(3,NFAC8)=velche(min(ivm1,LFAC(IAD+2)),min(ivm2,IEL)) 221 : xFAC8(4,NFAC8)=velche(min(ivm1,LFAC(IAD+3)),min(ivm2,IEL)) 222 : xFAC8(5,NFAC8)=velche(min(ivm1,LFAC(IAD+4)),min(ivm2,IEL)) 223 : xFAC8(6,NFAC8)=velche(min(ivm1,LFAC(IAD+5)),min(ivm2,IEL)) 224 : xFAC8(7,NFAC8)=velche(min(ivm1,LFAC(IAD+6)),min(ivm2,IEL)) 225 : xFAC8(8,NFAC8)=velche(min(ivm1,LFAC(IAD+7)),min(ivm2,IEL)) 226 : endif 227 : IFAC8(9,NFAC8)=0 228 : ipa=icpr(ifac8(1,nfac8)) 229 : ipb=icpr(ifac8(2,nfac8)) 230 : ipc=icpr(ifac8(3,nfac8)) 231 : ipd=icpr(ifac8(4,nfac8)) 232 : ipe=icpr(ifac8(5,nfac8)) 233 : ipf=icpr(ifac8(6,nfac8)) 234 : ipg=icpr(ifac8(7,nfac8)) 235 : iph=icpr(ifac8(8,nfac8)) 236 : IF (IVU(ipa).EQ.1.AND.IVU(ipb).EQ.1.AND.IVU(ipc).EQ.1 237 : > .AND.IVU(ipd).EQ.1.AND.IVU(ipe).EQ.1 238 : > .AND.IVU(ipf).EQ.1.AND.IVU(ipg).EQ.1 239 : > .AND.IVU(iph).EQ.1) 240 : > IFAC8(9,NFAC8)=1 241 : GOTO 67 242 : 66 NFACP=NFACP+1 243 : IFACP(15,NFACP)=1 244 : NSOMP(NFACP)= NBNN 245 : DO 660, IDC = 1, NBNN 246 : IFACP(IDC,NFACP)=IPT1.NUM(LFAC(IAD+IDC-1),IEL) 247 : if (mcham.ne.0) then 248 : XFACP(IDC,NFACP)=velche(min(ivm1,LFAC(IAD+IDC-1)), 249 : > min(ivm2,IEL)) 250 : endif 251 : IF (IVU(ICPR(IFACP(IDC,NFACP))).NE.1) IFACP(15,NFACP)=0 252 : 660 CONTINUE 253 : GOTO 67 254 : 67 CONTINUE 255 : 80 CONTINUE 256 : 60 CONTINUE 257 : IF (LISOUS(/1).NE.0) SEGDES IPT1 258 : goto 50 259 : 52 CONTINUE 260 : do 68 iel=1,nbelem 261 : if (ipt1.itypel.eq.2) then 262 : nseg2=nseg2+1 263 : ISEG2(1,NSEG2)=IPT1.NUM(1,IEL) 264 : ISEG2(2,NSEG2)=IPT1.NUM(2,IEL) 265 : elseif (ipt1.itypel.eq.3) then 266 : nseg3=nseg3+1 267 : ISEG3(1,NSEG3)=IPT1.NUM(1,IEL) 268 : ISEG3(2,NSEG3)=IPT1.NUM(2,IEL) 269 : ISEG3(3,NSEG3)=IPT1.NUM(3,IEL) 270 : endif 271 : if (mcham.ne.0) then 272 : if (ipt1.itypel.eq.2) then 273 : xSEG2(1,NSEG2)=velche(min(ivm1,1),min(ivm2,IEL)) 274 : xSEG2(2,NSEG2)=velche(min(ivm1,2),min(ivm2,IEL)) 275 : elseif (ipt1.itypel.eq.3) then 276 : xSEG3(1,NSEG3)=velche(min(ivm1,1),min(ivm2,IEL)) 277 : xSEG3(2,NSEG3)=velche(min(ivm1,2),min(ivm2,IEL)) 278 : xSEG3(3,NSEG3)=velche(min(ivm1,3),min(ivm2,IEL)) 279 : endif 280 : endif 281 : 68 continue 282 : IF (LISOUS(/1).NE.0) SEGDES IPT1 283 : 50 CONTINUE 284 : SEGDES MELEME 285 : C IF FAUT MAINTENANT RETASSER ET CLASSER LES TABLEAUX DES FACES 286 : C PROBLEME ELLES NE SONT PAS FORCEMENT DECRITES DE LA MEME FACON 287 : NFAC=NSEG2+NSEG3+NFAC3+NFAC4+NFAC6+NFAC8+NFACP 288 : SEGINI TFAC,KFAK 289 : IF (NFAC3.EQ.0) GOTO 101 290 : DO 102 I=1,NFAC3 291 : TFAC(I)=(XPROJ(3,ICPR(IFAC3(1,I)))+XPROJ(3,ICPR(IFAC3(2,I))) 292 : # + XPROJ(3,ICPR(IFAC3(3,I))))/3. 293 : IF (IFAC3(4,I).EQ.1) TFAC(I)=TFAC(I)-TDIST 294 : IF (IFAC3(4,I).EQ.2) TFAC(I)=TFAC(I)-2*TDIST 295 : KFAK(I)=I 296 : IF (IFAC3(4,I).EQ.3) KFAK(I)=0 297 : 102 CONTINUE 298 : 101 CONTINUE 299 : IF (NFAC4.EQ.0) GOTO 103 300 : DO 104 I=1,NFAC4 301 : ipa=icpr(ifac4(1,i)) 302 : ipb=icpr(ifac4(2,i)) 303 : ipc=icpr(ifac4(3,i)) 304 : ipd=icpr(ifac4(4,i)) 305 : TFAC(I+NFAC3)=(XPROJ(3,ipa)+xproj(3,ipb)+ 306 : > xproj(3,ipc)+xproj(3,ipd))/4.D0 307 : IF (IFAC4(5,I).EQ.1) TFAC(I+NFAC3)=TFAC(I+NFAC3)-TDIST 308 : KFAK(I+NFAC3)=I+NFAC 309 : 104 CONTINUE 310 : 103 CONTINUE 311 : IF (NFAC6.EQ.0) GOTO 105 312 : DO 106 I=1,NFAC6 313 : ipa=icpr(ifac6(1,i)) 314 : ipb=icpr(ifac6(2,i)) 315 : ipc=icpr(ifac6(3,i)) 316 : ipd=icpr(ifac6(4,i)) 317 : ipe=icpr(ifac6(5,i)) 318 : ipf=icpr(ifac6(6,i)) 319 : TFAC(I+NFAC3+NFAC4)=(XPROJ(3,ipa)+xproj(3,ipb)+ 320 : > xproj(3,ipc)+xproj(3,ipd)+xproj(3,ipe)+xproj(3,ipf) 321 : > )/6.D0 322 : IF (IFAC6(7,I).EQ.1) TFAC(I+NFAC3+NFAC4)=TFAC(I+NFAC3+NFAC4) 323 : # -TDIST 324 : KFAK(I+NFAC3+NFAC4)=I+2*NFAC 325 : 106 CONTINUE 326 : 105 CONTINUE 327 : IF (NFAC8.EQ.0) GOTO 107 328 : DO 108 I=1,NFAC8 329 : ipa=icpr(ifac8(1,i)) 330 : ipb=icpr(ifac8(2,i)) 331 : ipc=icpr(ifac8(3,i)) 332 : ipd=icpr(ifac8(4,i)) 333 : ipe=icpr(ifac8(5,i)) 334 : ipf=icpr(ifac8(6,i)) 335 : ipg=icpr(ifac8(7,i)) 336 : iph=icpr(ifac8(8,i)) 337 : TFAC(I+NFAC3+NFAC4+NFAC6)=(XPROJ(3,ipa)+xproj(3,ipb)+ 338 : > xproj(3,ipc)+xproj(3,ipd)+xproj(3,ipe)+xproj(3,ipf)+ 339 : > xproj(3,ipg)+xproj(3,iph))/8.D0 340 : IF (IFAC8(9,I).EQ.1) TFAC(I+NFAC3+NFAC4+NFAC6)= 341 : # TFAC(I+NFAC3+NFAC4+NFAC6)-TDIST 342 : KFAK(I+NFAC3+NFAC4+NFAC6)=I+3*NFAC 343 : 108 CONTINUE 344 : 107 CONTINUE 345 : IF (NSEG2.EQ.0) GOTO 109 346 : DO 110 I=1,NSEG2 347 : TFAC(I+NFAC3+NFAC4+NFAC6+NFAC8)=(XPROJ(3,ICPR(ISEG2(1,I)))+ 348 : # XPROJ(3,ICPR(ISEG2(2,I))))/2. 349 : KFAK(I+NFAC3+NFAC4+NFAC6+NFAC8)=I+4*NFAC 350 : 110 CONTINUE 351 : 109 CONTINUE 352 : IF (NSEG3.EQ.0) GOTO 111 353 : DO 112 I=1,NSEG3 354 : TFAC(I+NFAC3+NFAC4+NFAC6+NFAC8+NSEG2)=(XPROJ(3,ICPR(ISEG3(1,I)))+ 355 : # XPROJ(3,ICPR(ISEG3(2,I)))+XPROJ(3,ICPR(ISEG3(3,I))))/2. 356 : KFAK(I+NFAC3+NFAC4+NFAC6+NFAC8+NSEG2)=I+5*NFAC 357 : 112 CONTINUE 358 : 111 CONTINUE 359 : IF (NFACP.EQ.0) GOTO 121 360 : DO 122 I=1,NFACP 361 : XXXX = 0. 362 : DO 123 J=1, NSOMP(I) 363 : XXXX = XXXX + XPROJ(3,ICPR(IFACP(J,I))) 364 : 123 CONTINUE 365 : TFAC(I+NFAC3+NFAC4+NFAC6+NFAC8+NSEG2+NSEG3)= XXXX/NSOMP(I) 366 : KFAK(I+NFAC3+NFAC4+NFAC6+NFAC8+NSEG2+NSEG3)=I+6*NFAC 367 : 122 CONTINUE 368 : 121 CONTINUE 369 : C IL N'Y A PLUS QU'A TRIER ET RETASSER KFAK SUIVANT TFAC 370 : SEGINI NAUX 371 : IF (IDIM.EQ.2) GOTO 209 372 : NAUX(1)=1 373 : NAUX(2)=NFAC 374 : IZ=2 375 : 208 IZ=IZ-1 376 : IF (IZ.LE.0) GOTO 209 377 : IPB=NAUX(IZ*2-1) 378 : IPH=NAUX(IZ*2) 379 : IF(IPB.GE.IPH) GOTO 208 380 : JPB=IPB-1 381 : JPH=IPH+1 382 : C CALCUL DU PIVOT 383 : PV=0. 384 : * DO 207 J=IPB,IPH 385 : * PV=PV+TFAC(J) 386 : *207 CONTINUE 387 : * PV=PV/(IPH-IPB+1) 388 : PV=(TFAC(IPB)+TFAC(IPH))/2. 389 : 242 JPB=JPB+1 390 : IF (JPH.EQ.JPB) GOTO 245 391 : IF (TFAC(JPB).LT.PV) GOTO 243 392 : GOTO 242 393 : 243 JPH=JPH-1 394 : IF (JPH.EQ.JPB) GOTO 245 395 : IF (TFAC(JPH).GT.PV) GOTO 244 396 : GOTO 243 397 : 244 IAUX=KFAK(JPB) 398 : KFAK(JPB)=KFAK(JPH) 399 : KFAK(JPH)=IAUX 400 : TAUX=TFAC(JPB) 401 : TFAC(JPB)=TFAC(JPH) 402 : TFAC(JPH)=TAUX 403 : GOTO 242 404 : 245 IF (JPB.GE.IPB) GOTO 247 405 : JPB=JPB+1 406 : JPH=JPH+2 407 : GOTO 248 408 : 247 IF (JPH.LE.IPH) GOTO 249 409 : JPB=JPB-2 410 : JPH=JPH-1 411 : GOTO 248 412 : 249 IF (TFAC(JPB).LE.PV) GOTO 250 413 : IF (JPH.EQ.IPH) GOTO 251 414 : 252 JPH=JPH+1 415 : GOTO 248 416 : 250 IF (JPB.EQ.IPB) GOTO 252 417 : 251 JPB=JPB-1 418 : 248 IF (JPB.EQ.IPB) GOTO 253 419 : NAUX(2*IZ)=JPB 420 : IZ=IZ+1 421 : 253 IF (JPH.EQ.IPH) GOTO 208 422 : NAUX(2*IZ)=IPH 423 : NAUX(2*IZ-1)=JPH 424 : IZ=IZ+1 425 : GOTO 208 426 : 209 CONTINUE 427 : C LES FACES SONT CLASSEES DANS KFAK IL AURAIT FALLU ELIMINER LES FACES 428 : C EN DOUBLE . IL NE RESTE PLUS QU'A TRACER 429 : DO 300 I=1,NFAC 430 : IFF=KFAK(I) 431 : IF (IFF.EQ.0) GOTO 300 432 : IT=(IFF-1)/NFAC 433 : IF=IFF-IT*NFAC 434 : IT=IT+1 435 : IOK=0 436 : GOTO (310,320,330,340,350,360,370),IT 437 : 310 CONTINUE 438 : IOK=0 439 : IPA=ICPR(IFAC3(1,IF)) 440 : IF (IVU(IPA).EQ.1) IOK=1 441 : XX(1)=XPROJ(1,IPA) 442 : YY(1)=XPROJ(2,IPA) 443 : ZZ(1)=XPROJ(3,IPA) 444 : if (mcham.eq.0) then 445 : VV(1)=VCPCHA(IFAC3(1,IF)) 446 : else 447 : VV(1)=xFAC3(1,IF) 448 : endif 449 : IPB=ICPR(IFAC3(2,IF)) 450 : IF (IVU(IPB).EQ.1) IOK=1 451 : XX(2)=XPROJ(1,IPB) 452 : YY(2)=XPROJ(2,IPB) 453 : ZZ(2)=XPROJ(3,IPB) 454 : if (mcham.eq.0) then 455 : VV(2)=VCPCHA(IFAC3(2,IF)) 456 : else 457 : VV(2)=xFAC3(2,IF) 458 : endif 459 : IPC=ICPR(IFAC3(3,IF)) 460 : IF (IVU(IPC).EQ.1) IOK=1 461 : XX(3)=XPROJ(1,IPC) 462 : YY(3)=XPROJ(2,IPC) 463 : ZZ(3)=XPROJ(3,IPC) 464 : if (mcham.eq.0) then 465 : VV(3)=VCPCHA(IFAC3(3,IF)) 466 : else 467 : VV(3)=xFAC3(3,IF) 468 : endif 468 : IF (IOK.EQ.1) CALL TRISO(VCHC,XX,YY,ZZ,VV,3,NISO) 469 : GOTO 380 470 : 320 continue 471 : * cas des quadrangles on rajoute un point au milieu 472 : ia=icpr(ifac4(1,if)) 473 : ib=icpr(ifac4(2,if)) 474 : ic=icpr(ifac4(3,if)) 475 : id=icpr(ifac4(4,if)) 476 : IF (IVU(IA).NE.1.AND.IVU(IB).NE.1.AND.IVU(IC).NE.1.AND. 477 : * IVU(Id).NE.1) goto 380 478 : xxm=(XPROJ(1,IA)+XPROJ(1,IB)+XPROJ(1,IC)+XPROJ(1,ID))/4. 479 : yym=(XPROJ(2,IA)+XPROJ(2,IB)+XPROJ(2,IC)+XPROJ(2,ID))/4. 480 : zzm=(XPROJ(3,IA)+XPROJ(3,IB)+XPROJ(3,IC)+XPROJ(3,ID))/4. 481 : if (mcham.eq.0) then 482 : vvm=(VCPCHA(ifac4(1,if))+VCPCHA(ifac4(2,if))+ 483 : * VCPCHA(ifac4(3,if))+VCPCHA(ifac4(4,if)))/4. 484 : else 485 : vvm=(xfac4(1,if)+xfac4(2,if)+ 486 : * xfac4(3,if)+xfac4(4,if))/4. 487 : endif 488 : xx(1)=xxm 489 : yy(1)=yym 490 : zz(1)=zzm 491 : VV(1)=vvm 492 : xx(2)=xproj(1,ia) 493 : yy(2)=xproj(2,ia) 494 : zz(2)=xproj(3,ia) 495 : xx(3)=xproj(1,ib) 496 : yy(3)=xproj(2,ib) 497 : zz(3)=xproj(3,ib) 498 : xx(4)=xproj(1,ic) 499 : yy(4)=xproj(2,ic) 500 : zz(4)=xproj(3,ic) 501 : xx(5)=xproj(1,id) 502 : yy(5)=xproj(2,id) 503 : zz(5)=xproj(3,id) 504 : if (mcham.eq.0) then 505 : VV(2)=VCPCHA(ifac4(1,if)) 506 : VV(3)=VCPCHA(ifac4(2,if)) 507 : VV(4)=VCPCHA(ifac4(3,if)) 508 : VV(5)=VCPCHA(ifac4(4,if)) 509 : else 510 : VV(2)=xfac4(1,if) 511 : VV(3)=xfac4(2,if) 512 : VV(4)=xfac4(3,if) 513 : VV(5)=xfac4(4,if) 514 : endif 514 : CALL TCISO(VCHC,XX,YY,ZZ,VV,5,NISO) 515 : goto 380 516 : 517 : 330 continue 518 : * cas des triangles a 6 noeuds on 519 : ipa=icpr(ifac6(1,if)) 520 : ipb=icpr(ifac6(2,if)) 521 : ipc=icpr(ifac6(3,if)) 522 : ipd=icpr(ifac6(4,if)) 523 : ipe=icpr(ifac6(5,if)) 524 : ipf=icpr(ifac6(6,if)) 525 : IF (IVU(IpA).NE.1.AND.IVU(IpB).NE.1.AND.IVU(IpC).NE.1.AND. 526 : * IVU(IpD).NE.1.AND.IVU(IpE).NE.1.AND.IVU(IpF).NE.1) goto 380 527 : * cas du quadrangle on rajoute un point central 528 : xxm=(XPROJ(1,IpB)+XPROJ(1,IpD)+XPROJ(1,IpF))/3 529 : yym=(XPROJ(2,IpB)+XPROJ(2,IpD)+XPROJ(2,IpF))/3 530 : zzm=(XPROJ(3,IpB)+XPROJ(3,IpD)+XPROJ(3,IpF))/3 531 : if (mcham.eq.0) then 532 : vvm=(VCPCHA(ifac6(2,if))+VCPCHA(ifac6(4,if))+ 533 : > VCPCHA(ifac6(6,if)))/3 534 : else 535 : vvm=(xfac6(2,if)+xfac6(4,if)+xfac6(6,if))/3 536 : endif 537 : * 538 : xx(1)=xxm 539 : yy(1)=yym 540 : zz(1)=zzm 541 : VV(1)=vvm 542 : xx(2)=xproj(1,ipa) 543 : yy(2)=xproj(2,ipa) 544 : zz(2)=xproj(3,ipa) 545 : xx(3)=xproj(1,ipb) 546 : yy(3)=xproj(2,ipb) 547 : zz(3)=xproj(3,ipb) 548 : xx(4)=xproj(1,ipc) 549 : yy(4)=xproj(2,ipc) 550 : zz(4)=xproj(3,ipc) 551 : xx(5)=xproj(1,ipd) 552 : yy(5)=xproj(2,ipd) 553 : zz(5)=xproj(3,ipd) 554 : xx(6)=xproj(1,ipe) 555 : yy(6)=xproj(2,ipe) 556 : zz(6)=xproj(3,ipe) 557 : xx(7)=xproj(1,ipf) 558 : yy(7)=xproj(2,ipf) 559 : zz(7)=xproj(3,ipf) 560 : if (mcham.eq.0) then 561 : VV(2)=VCPCHA(ifac6(1,if)) 562 : VV(3)=VCPCHA(ifac6(2,if)) 563 : VV(4)=VCPCHA(ifac6(3,if)) 564 : VV(5)=VCPCHA(ifac6(4,if)) 565 : VV(6)=VCPCHA(ifac6(5,if)) 566 : VV(7)=VCPCHA(ifac6(6,if)) 567 : else 568 : VV(2)=xfac6(1,if) 569 : VV(3)=xfac6(2,if) 570 : VV(4)=xfac6(3,if) 571 : VV(5)=xfac6(4,if) 572 : VV(6)=xfac6(5,if) 573 : VV(7)=xfac6(6,if) 574 : endif 574 : CALL TCISO(VCHC,XX,YY,ZZ,VV,7,NISO) 575 : goto 380 576 : 340 continue 577 : * cas des quadrangles a 8 noeuds 578 : ipa=icpr(ifac8(1,if)) 579 : ipb=icpr(ifac8(2,if)) 580 : ipc=icpr(ifac8(3,if)) 581 : ipd=icpr(ifac8(4,if)) 582 : ipe=icpr(ifac8(5,if)) 583 : ipf=icpr(ifac8(6,if)) 584 : ipg=icpr(ifac8(7,if)) 585 : iph=icpr(ifac8(8,if)) 586 : IF (IVU(IpA).NE.1.AND.IVU(IpB).NE.1.AND.IVU(IpC).NE.1.AND. 587 : * IVU(IpD).NE.1.AND.IVU(IpE).NE.1.AND.IVU(IpF).NE.1.AND. 588 : > IVU(IpG).NE.1.AND.IVU(IpH).NE.1) goto 380 589 : * cas du quadrangle on rajoute un point central 590 : xxm=-0.25*XPROJ(1,IpA)+0.5*XPROJ(1,IpB) 591 : * -0.25*XPROJ(1,IpC)+0.5*XPROJ(1,IpD) 592 : * -0.25*XPROJ(1,IpE)+0.5*XPROJ(1,IpF) 593 : * -0.25*XPROJ(1,IpG)+0.5*XPROJ(1,IpH) 594 : yym=-0.25*XPROJ(2,IpA)+0.5*XPROJ(2,IpB) 595 : * -0.25*XPROJ(2,IpC)+0.5*XPROJ(2,IpD) 596 : * -0.25*XPROJ(2,IpE)+0.5*XPROJ(2,IpF) 597 : * -0.25*XPROJ(2,IpG)+0.5*XPROJ(2,IpH) 598 : zzm=-0.25*XPROJ(3,IpA)+0.5*XPROJ(3,IpB) 599 : * -0.25*XPROJ(3,IpC)+0.5*XPROJ(3,IpD) 600 : * -0.25*XPROJ(3,IpE)+0.5*XPROJ(3,IpF) 601 : * -0.25*XPROJ(3,IpG)+0.5*XPROJ(3,IpH) 602 : if (mcham.eq.0) then 603 : vvm= 604 : * -0.25*VCPCHA(ifac8(1,if))+0.5*VCPCHA(ifac8(2,if)) 605 : * -0.25*VCPCHA(ifac8(3,if))+0.5*VCPCHA(ifac8(4,if)) 606 : vvm = vvm 607 : * -0.25*VCPCHA(ifac8(5,if))+0.5*VCPCHA(ifac8(6,if)) 608 : * -0.25*VCPCHA(ifac8(7,if))+0.5*VCPCHA(ifac8(8,if)) 609 : else 610 : vvm= 611 : * -0.25*xfac8(1,if)+0.5*xfac8(2,if) 612 : * -0.25*xfac8(3,if)+0.5*xfac8(4,if) 613 : * -0.25*xfac8(5,if)+0.5*xfac8(6,if) 614 : * -0.25*xfac8(7,if)+0.5*xfac8(8,if) 615 : endif 616 : * 617 : xx(1)=xxm 618 : yy(1)=yym 619 : zz(1)=zzm 620 : VV(1)=vvm 621 : xx(2)=xproj(1,ipa) 622 : yy(2)=xproj(2,ipa) 623 : zz(2)=xproj(3,ipa) 624 : xx(3)=xproj(1,ipb) 625 : yy(3)=xproj(2,ipb) 626 : zz(3)=xproj(3,ipb) 627 : xx(4)=xproj(1,ipc) 628 : yy(4)=xproj(2,ipc) 629 : zz(4)=xproj(3,ipc) 630 : xx(5)=xproj(1,ipd) 631 : yy(5)=xproj(2,ipd) 632 : zz(5)=xproj(3,ipd) 633 : xx(6)=xproj(1,ipe) 634 : yy(6)=xproj(2,ipe) 635 : zz(6)=xproj(3,ipe) 636 : xx(7)=xproj(1,ipf) 637 : yy(7)=xproj(2,ipf) 638 : zz(7)=xproj(3,ipf) 639 : xx(8)=xproj(1,ipg) 640 : yy(8)=xproj(2,ipg) 641 : zz(8)=xproj(3,ipg) 642 : xx(9)=xproj(1,iph) 643 : yy(9)=xproj(2,iph) 644 : zz(9)=xproj(3,iph) 645 : if (mcham.eq.0) then 646 : VV(2)=VCPCHA(ifac8(1,if)) 647 : VV(3)=VCPCHA(ifac8(2,if)) 648 : VV(4)=VCPCHA(ifac8(3,if)) 649 : VV(5)=VCPCHA(ifac8(4,if)) 650 : VV(6)=VCPCHA(ifac8(5,if)) 651 : VV(7)=VCPCHA(ifac8(6,if)) 652 : VV(8)=VCPCHA(ifac8(7,if)) 653 : VV(9)=VCPCHA(ifac8(8,if)) 654 : else 655 : VV(2)=xfac8(1,if) 656 : VV(3)=xfac8(2,if) 657 : VV(4)=xfac8(3,if) 658 : VV(5)=xfac8(4,if) 659 : VV(6)=xfac8(5,if) 660 : VV(7)=xfac8(6,if) 661 : VV(8)=xfac8(7,if) 662 : VV(9)=xfac8(8,if) 663 : endif 663 : CALL TCISO(VCHC,XX,YY,ZZ,VV,9,NISO) 664 : goto 380 665 : 350 CONTINUE 666 : * cas du seg2 667 : IOK=0 668 : IPA=ICPR(ISEG2(1,IF)) 669 : IF (IVU(IPA).EQ.1) IOK=1 670 : XX(1)=XPROJ(1,IPA) 671 : YY(1)=XPROJ(2,IPA) 672 : ZZ(1)=XPROJ(3,IPA) 673 : if (mcham.eq.0) then 674 : VV(1)=VCPCHA(ISEG2(1,IF)) 675 : else 676 : VV(1)=xseg2(1,IF) 677 : endif 678 : IPB=ICPR(Iseg2(2,IF)) 679 : IF (IVU(IPB).EQ.1) IOK=1 680 : XX(2)=XPROJ(1,IPB) 681 : YY(2)=XPROJ(2,IPB) 682 : ZZ(2)=XPROJ(3,IPB) 683 : if (mcham.eq.0) then 684 : VV(2)=VCPCHA(Iseg2(2,IF)) 685 : else 686 : VV(2)=xseg2(2,IF) 687 : endif 687 : IF (IOK.EQ.1) CALL TRISO(VCHC,XX,YY,ZZ,VV,2,NISO) 688 : GOTO 380 689 : 360 CONTINUE 690 : * cas du seg3 691 : IOK=0 692 : IPA=ICPR(ISEG3(1,IF)) 693 : IF (IVU(IPA).EQ.1) IOK=1 694 : XX(1)=XPROJ(1,IPA) 695 : YY(1)=XPROJ(2,IPA) 696 : ZZ(1)=XPROJ(3,IPA) 697 : if (mcham.eq.0) then 698 : VV(1)=VCPCHA(ISEG3(1,IF)) 699 : else 700 : VV(1)=xseg3(1,IF) 701 : endif 702 : IPB=ICPR(Iseg3(2,IF)) 703 : IF (IVU(IPB).EQ.1) IOK=1 704 : XX(2)=XPROJ(1,IPB) 705 : YY(2)=XPROJ(2,IPB) 706 : ZZ(2)=XPROJ(3,IPB) 707 : if (mcham.eq.0) then 708 : VV(2)=VCPCHA(Iseg3(2,IF)) 709 : else 710 : VV(2)=xseg3(2,IF) 711 : endif 711 : IF (IOK.EQ.1) CALL TRISO(VCHC,XX,YY,ZZ,VV,2,NISO) 712 : IOK=0 713 : IPA=ICPR(ISEG3(2,IF)) 714 : IF (IVU(IPA).EQ.1) IOK=1 715 : XX(1)=XPROJ(1,IPA) 716 : YY(1)=XPROJ(2,IPA) 717 : ZZ(1)=XPROJ(3,IPA) 718 : if (mcham.eq.0) then 719 : VV(1)=VCPCHA(ISEG3(2,IF)) 720 : else 721 : VV(1)=xseg3(2,IF) 722 : endif 723 : IPB=ICPR(Iseg3(3,IF)) 724 : IF (IVU(IPB).EQ.1) IOK=1 725 : XX(2)=XPROJ(1,IPB) 726 : YY(2)=XPROJ(2,IPB) 727 : ZZ(2)=XPROJ(3,IPB) 728 : if (mcham.eq.0) then 729 : VV(2)=VCPCHA(Iseg3(3,IF)) 730 : else 731 : VV(2)=xseg3(3,IF) 732 : endif 732 : IF (IOK.EQ.1) CALL TRISO(VCHC,XX,YY,ZZ,VV,2,NISO) 733 : GOTO 380 734 : 370 CONTINUE 735 : * 736 : * POLYGONE 737 : * 738 : IOK=1 739 : NBCOT = NSOMP(IF) 740 : * 741 : * Recherche des coordonnees et valeurs au centre du polygone 742 : * 743 : XXM = 0. 744 : YYM = 0. 745 : VVM = 0. 746 : IVUE = 1 747 : DO 3700 ICT = 1, NBCOT 748 : * 749 : NUPT = IFACP(ICT, IF) 750 : IDPT = ICPR(NUPT) 751 : XXM = XPROJ(1,IDPT) + XXM 752 : YYM = XPROJ(2,IDPT) + YYM 753 : IF (MCHAM.EQ.0) THEN 754 : VVM = VCPCHA(NUPT) + VVM 755 : ELSE 756 : VVM = XFACP(ICT, IF) + VVM 757 : ENDIF 758 : IF (IVU(IDPT).NE.1) IVUE = 0 759 : * 760 : 3700 CONTINUE 761 : * 762 : IF (IVUE.EQ.1) THEN 763 : * 764 : XX(3)=XXM/NBCOT 765 : YY(3)=YYM/NBCOT 766 : VV(3)=VVM/NBCOT 767 : * 768 : NUPT2 = IFACP(NBCOT, IF) 769 : ID2 = ICPR(NUPT2) 770 : XX(2)=XPROJ(1,ID2) 771 : YY(2)=XPROJ(2,ID2) 772 : ZZ(2)=XPROJ(3,ID2) 773 : IF (MCHAM.EQ.0) THEN 774 : VV(2) = VCPCHA(NUPT2) 775 : ELSE 776 : VV(2) = XFACP(NBCOT, IF) 777 : ENDIF 778 : * 779 : * Boucle sur tous les triangles 780 : * 781 : DO 3705 , ICT = 1, NBCOT 782 : * 783 : NUPT1 = IFACP(ICT, IF) 784 : ID1 = ICPR(NUPT1) 785 : * 786 : XX(1)=XPROJ(1,ID1) 787 : YY(1)=XPROJ(2,ID1) 788 : ZZ(1)=XPROJ(3,ID1) 789 : IF (MCHAM.EQ.0) THEN 790 : VV(1)= VCPCHA(NUPT1) 791 : ELSE 792 : VV(1) = XFACP(ICT, IF) 793 : ENDIF 793 : CALL TRISO(VCHC,XX,YY,ZZ,VV,3,NISO) 794 : XX(2)= XX(1) 795 : YY(2)= YY(1) 796 : ZZ(2)= ZZ(1) 797 : VV(2)= VV(1) 798 : * 799 : 3705 CONTINUE 800 : * 801 : ENDIF 802 : * 803 : GOTO 380 804 : 380 CONTINUE 805 : 300 CONTINUE 806 : C 'EST FINI 807 : SEGACT MELEME 808 : IF (LISOUS(/1).NE.0) THEN 809 : NBSOUS=LISOUS(/1) 810 : * IF (MCOUP.NE.0) NBSOUS=NBSOUS-1 811 : DO 490 IO=1,NBSOUS 812 : IPT2=LISOUS(IO) 813 : segact ipt2 814 : if (ipt2.itypel.gt.3.AND.ipt2.itypel.NE.32) then 815 : SEGSUP IPT2 816 : if (mcham.ne.0) then 817 : melval=lisref(io) 818 : segsup melval 819 : endif 820 : else 821 : segdes ipt2 822 : endif 823 : 490 CONTINUE 824 : ENDIF 825 : if (itypel.eq.0) SEGSUP MELEME 826 : MELEME=MELSAU 827 : SEGSUP ISEG2,ISEG3,IFAC3,IFAC4,IFAC6,IFAC8,NAUX,TFAC,KFAK 828 : if (mcham.ne.0) SEGSUP XSEG2,XSEG3,xFAC3,xFAC4,xFAC6,xFAC8 829 : RETURN 830 : END 831 : 832 : 833 : 834 : 835 : 836 : 837 : 838 : 839 : 840 : 841 : 842 : 843 : 844 : 845 : 846 : 847 : 848 : 849 :
© Cast3M 2003 - All rights reserved.
Disclaimer