Note a:
This was formerly the 3D printer problem. My
apologies if this was already solved. This code is
(unproven) a general solution for tetrahedral
meshes.
)
Note 'MESH and algorithm'
A mesh is a boxed vector length 2.
0{::MESH is a rank 2 array of indexes into the nodes.
1{::MESH is a rank 2 array of nodes.
Each row of indexes defines an element.
Each row of nodes contains coordinates of a node.
Consider simplexes in 3 dimensions; tetrahedra.
Number them in a right-handed fashion.
Wrap your right hand around the first three nodes,
your thumb points into the half space that includes
the 4th vertex.
Example with a single tetrahedron:
]MESH=: (i. 1 4) ; 0 , (, 1) , 1 1 ,: 1 1 1
+-------+-----+
|0 1 2 3|0 0 0|
| |1 0 0|
| |1 1 0|
| |1 1 1|
+-------+-----+
I propose that the parallelepiped volumes formed by
each face and an additional point are positive iff
that additional point is in the tetrahedron. Proof?
This code hinges upon this missing proof.
The right-handed face indexes into the tetrahedron
vertices are TETRAHEDRON_FACES
)
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
