I think your theorem is true.
Same idea for in_tetrahedra, but using determinants to check sign:
int=:4 :'(-: 1&|.)@(0<!.0 *&_1 1 _1 1)@(1&((-/ .*)\.))"2 ({ -"1&x)&>/ y'
Use intolerant comparison for sign checks to help make unique assignment
to cells on boundaries.
Henry Rich
On 7/3/2012 10:02 PM, David Ward Lambert wrote:
> F=: TETRAHEDRON_FACES=: ".;._2(0 :0)
> 0 3 1
> 1 3 2
> 2 3 0
> 0 1 2
> )
>
> NB. construct a test mesh
> E=: 0 1 2 3 ; 5 1 2 3 ; 0 2 1 3 ; 0 0 2 3 ; 0 1 2 4
> N=: 0 0 0 ; 1 0 0 ; 1 1 0 ; 1 1 1 ; 0 1 0 ; 1r2 0 0
> MESH=: E ;&:> N NB. only the first 2 elements are valid.
>
> mp=: $:~ : (+/ .*)
> assert 25 -: mp 3 4
> assert 0 -: 1 0 0 mp 0 1 0
>
> crossproduct=: 1 |. ([ * 1 |. ]) - ] * 1 |. [ NB. j forum
> assert 0 0 1 -: 1 0 0 crossproduct 0 1 0
> assert 0 0 _1 -: 1 0 0 crossproduct~ 0 1 0
>
> NB. monad parallelepiped_volume
> NB. 4 3 -: $ y
> ppv=: ({: mp [: crossproduct/ 2&{.)@(2&(-~/\))"2 :[:
> parallelepiped_volume=: ppv
> assert 1 -: ppv 4 {. 1 {:: MESH NB. 1st 4 nodes -: unit cube
>
> Note 'verify_tetrahedra'
> NB. y a MESH, return Boolean vector of the valid tetrahedra
> verify_tetrahedra=: 3 : 0
> 'E N'=. y
> 0 < parallelepiped_volume E { N
> )
> verify_tetrahedra=: 0 < parallelepiped_volume@:{&:>/
> assert 1 1 0 0 0 -: verify_tetrahedra MESH
>
> NB. in_tetrahedra returns Boolean vector of tets containing x
> NB. x is a point, y a mesh
> NB. Discard bad elements sooner in a production application.
> in_tetrahedra=: dyad define
> 'E N'=. y
> BASES=. TETRAHEDRON_FACES {"_ 1 E
> VERTEXES=. x ,"2~ BASES { N
> (verify_tetrahedra y) *. 0 (*./ .<:)"1 ppv VERTEXES
> )
> assert 1 0 0 0 0 -: 0.4 0.4 0.1 in_tetrahedra MESH
> assert 1 1 0 0 0 -: 0.7 0.1 0.1 in_tetrahedra MESH
> assert 0 0 0 0 0 -: (-1 1 1) in_tetrahedra MESH
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
