C VERVOL SOURCE JC220346 16/11/29 21:15:40 9221 C verification que le pt ip n'est pas dans un volume C ayant une arete ipi ipj C LOGICAL FUNCTION VERVOL(IPT,IP1,IP2,IP3,IP4,IP5,IP6,IP7,IP8) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC TDEMAIT dimension ipp(8),ivtemp(50) np=0 if (ip1.gt.0) then np=np+1 ipp(np)=ip1 endif if (ip2.gt.0) then np=np+1 ipp(np)=ip2 endif if (ip3.gt.0) then np=np+1 ipp(np)=ip3 endif if (ip4.gt.0) then np=np+1 ipp(np)=ip4 endif if (ip5.gt.0) then np=np+1 ipp(np)=ip5 endif if (ip6.gt.0) then np=np+1 ipp(np)=ip6 endif if (ip7.gt.0) then np=np+1 ipp(np)=ip7 endif if (ip8.gt.0) then np=np+1 ipp(np)=ip8 endif vervol=.true. * return * * boucle sur les elements * nvt=0 do 9 iq=1,np ip=ipp(iq) do 9 ipf=1,40 if=npf(ipf,iq) if (if.eq.0) goto 5 do 6 ifv=1,2 iv=nfv(ifv,if) if (iv.eq.0) goto 9 do 8 ivt=1,nvt if (ivtemp(ivt).eq.iv) goto 7 8 continue nvt=nvt+1 ivtemp(nvt)=iv 7 continue 6 continue 5 continue 9 continue do 10 ivt=1,nvt iv=ivtemp(ivt) it=0 do 20 ip=1,8 if (ivol(ip,iv).eq.ip1) it=it+1 if (ivol(ip,iv).eq.ip2) it=it+1 if (ivol(ip,iv).eq.ip3) it=it+1 if (ivol(ip,iv).eq.ip4) it=it+1 if (ivol(ip,iv).eq.ip5) it=it+1 if (ivol(ip,iv).eq.ip6) it=it+1 if (ivol(ip,iv).eq.ip7) it=it+1 if (ivol(ip,iv).eq.ip8) it=it+1 20 continue if (it.lt.2) goto 10 * l'element a 2 pt commun avec le notre ==> test supplementaire * cas du tetraedre if (ivol(9,iv).eq.25) then v1=vol(ipt,ivol(1,iv),ivol(2,iv),ivol(3,iv)) v2=vol(ipt,ivol(2,iv),ivol(1,iv),ivol(4,iv)) v3=vol(ipt,ivol(3,iv),ivol(2,iv),ivol(4,iv)) v4=vol(ipt,ivol(1,iv),ivol(3,iv),ivol(4,iv)) vv=v1+v2+v3+v4 * write (6,*) ' vervol vv v1 v2 v3 v4 ',vv,v1,v2,v3,v4 if (vv*v1.gt.0..and.vv*v2.gt.0..and.vv*v3.gt.0..and. * vv*v4.gt.0.) vervol=.false. endif if (.not.vervol) then * write (6,*) ' vervol element incorrect tetraedre' * write (6,*) ' vervol point teste ',ipt * write (6,*) xyz(1,ipt),xyz(2,ipt),xyz(3,ipt) ipv=ivol(1,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(2,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(3,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(4,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) return endif * cas de la pyramide if (ivol(9,iv).eq.35) then v1=vol(ipt,ivol(1,iv),ivol(2,iv),ivol(3,iv)) v2=vol(ipt,ivol(1,iv),ivol(3,iv),ivol(4,iv)) v3=vol(ipt,ivol(1,iv),ivol(2,iv),ivol(4,iv)) v4=vol(ipt,ivol(2,iv),ivol(3,iv),ivol(4,iv)) v5=vol(ipt,ivol(2,iv),ivol(1,iv),ivol(5,iv)) v6=vol(ipt,ivol(3,iv),ivol(2,iv),ivol(5,iv)) v7=vol(ipt,ivol(4,iv),ivol(3,iv),ivol(5,iv)) v8=vol(ipt,ivol(1,iv),ivol(4,iv),ivol(5,iv)) vv=v1+v2+v3+v4+v5+v6+v7+v8 if (vv*v1.gt.0..and.vv*v2.gt.0..and.vv*v3.gt.0..and. * vv*v4.gt.0..and.vv*v5.gt.0..and.vv*v6.gt.0..and. * vv*v7.gt.0..and.vv*v8.gt.0.) vervol=.false. endif if (.not.vervol) then * write (6,*) ' vervol element incorrect pyramide' * write (6,*) xyz(1,ipt),xyz(2,ipt),xyz(3,ipt) ipv=ivol(1,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(2,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(3,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(4,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(5,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) return endif * cas du prisme if (ivol(9,iv).eq.30) then v1=vol(ipt,ivol(1,iv),ivol(2,iv),ivol(3,iv)) v2=vol(ipt,ivol(1,iv),ivol(4,iv),ivol(5,iv)) v3=vol(ipt,ivol(1,iv),ivol(5,iv),ivol(2,iv)) v4=vol(ipt,ivol(1,iv),ivol(4,iv),ivol(2,iv)) v5=vol(ipt,ivol(4,iv),ivol(5,iv),ivol(2,iv)) v6=vol(ipt,ivol(2,iv),ivol(5,iv),ivol(6,iv)) v7=vol(ipt,ivol(2,iv),ivol(6,iv),ivol(3,iv)) v8=vol(ipt,ivol(2,iv),ivol(5,iv),ivol(3,iv)) v9=vol(ipt,ivol(5,iv),ivol(6,iv),ivol(3,iv)) v10=vol(ipt,ivol(3,iv),ivol(6,iv),ivol(4,iv)) v11=vol(ipt,ivol(3,iv),ivol(4,iv),ivol(1,iv)) v12=vol(ipt,ivol(3,iv),ivol(6,iv),ivol(1,iv)) v13=vol(ipt,ivol(4,iv),ivol(1,iv),ivol(2,iv)) v14=vol(ipt,ivol(6,iv),ivol(5,iv),ivol(4,iv)) vv=v1+v2+v3+v4+v5+v6+v7+v8+v9+v10+v11+v12+v13+v14 if (vv*v1.gt.0..and.vv*v2.gt.0..and.vv*v3.gt.0..and. * vv*v4.gt.0..and.vv*v5.gt.0..and.vv*v6.gt.0..and. * vv*v7.gt.0..and.vv*v8.gt.0..and.vv*v9.gt.0..and. * vv*v10.gt.0..and.vv*v11.gt.0..and.vv*v12.gt.0..and. * vv*v13.gt.0..and.vv*v14.gt.0.) vervol=.false. endif if (.not.vervol) then * write (6,*) ' vervol element incorrect prisme ' * write (6,*) ' vervol point teste ',ipt IF (IVERB.EQ.1) write (6,*) xyz(1,ipt),xyz(2,ipt),xyz(3,ipt) ipv=ivol(1,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(2,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(3,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(4,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(5,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(6,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) return endif * cas du cube if (ivol(9,iv).eq.20) then v1=vol(ipt,ivol(1,iv),ivol(2,iv),ivol(3,iv)) v2=vol(ipt,ivol(1,iv),ivol(3,iv),ivol(4,iv)) v3=vol(ipt,ivol(1,iv),ivol(2,iv),ivol(4,iv)) v4=vol(ipt,ivol(2,iv),ivol(3,iv),ivol(4,iv)) v5=vol(ipt,ivol(1,iv),ivol(5,iv),ivol(6,iv)) v6=vol(ipt,ivol(1,iv),ivol(6,iv),ivol(2,iv)) v7=vol(ipt,ivol(1,iv),ivol(5,iv),ivol(2,iv)) v8=vol(ipt,ivol(5,iv),ivol(6,iv),ivol(2,iv)) v9=vol(ipt,ivol(2,iv),ivol(6,iv),ivol(7,iv)) v10=vol(ipt,ivol(2,iv),ivol(7,iv),ivol(3,iv)) v11=vol(ipt,ivol(2,iv),ivol(6,iv),ivol(3,iv)) v12=vol(ipt,ivol(6,iv),ivol(7,iv),ivol(3,iv)) v13=vol(ipt,ivol(3,iv),ivol(7,iv),ivol(8,iv)) v14=vol(ipt,ivol(3,iv),ivol(8,iv),ivol(4,iv)) v15=vol(ipt,ivol(3,iv),ivol(7,iv),ivol(4,iv)) v16=vol(ipt,ivol(7,iv),ivol(8,iv),ivol(4,iv)) v17=vol(ipt,ivol(4,iv),ivol(8,iv),ivol(5,iv)) v18=vol(ipt,ivol(4,iv),ivol(5,iv),ivol(1,iv)) v19=vol(ipt,ivol(4,iv),ivol(8,iv),ivol(1,iv)) v20=vol(ipt,ivol(8,iv),ivol(5,iv),ivol(1,iv)) v21=vol(ipt,ivol(5,iv),ivol(8,iv),ivol(7,iv)) v22=vol(ipt,ivol(5,iv),ivol(7,iv),ivol(6,iv)) v23=vol(ipt,ivol(5,iv),ivol(8,iv),ivol(6,iv)) v24=vol(ipt,ivol(8,iv),ivol(7,iv),ivol(6,iv)) vv=v1+v2+v3+v4+v5+v6+v7+v8+v9+v10+v11+v12+ * v13+v14+v15+v16+v17+v18+v19+v20+v21+v22+v23+v24 if (vv*v1.gt.0..and.vv*v2.gt.0..and.vv*v3.gt.0..and. * vv*v4.gt.0..and.vv*v5.gt.0..and.vv*v6.gt.0..and. * vv*v7.gt.0..and.vv*v8.gt.0..and.vv*v9.gt.0..and. * vv*v10.gt.0..and.vv*v11.gt.0..and.vv*v12.gt.0..and. * vv*v13.gt.0..and.vv*v14.gt.0..and.vv*v15.gt.0..and. * vv*v16.gt.0..and.vv*v17.gt.0..and.vv*v18.gt.0..and. * vv*v19.gt.0..and.vv*v20.gt.0..and.vv*v21.gt.0..and. * vv*v22.gt.0..and.vv*v23.gt.0..and.vv*v24.gt.0.) * vervol=.false. endif if (.not.vervol) then * write (6,*) ' vervol element incorrect cube' * write (6,*) ' vervol point teste ',ipt * write (6,*) xyz(1,ipt),xyz(2,ipt),xyz(3,ipt) ipv=ivol(1,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(2,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(3,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(4,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(5,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(6,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(7,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) ipv=ivol(8,iv) * write (6,*) ' pt du volume ',ipv * write (6,*) xyz(1,ipv),xyz(2,ipv),xyz(3,ipv) return endif 10 continue end