I have no solution at present, but I think that turning to barycentric coordinates (or tetrahedral coordinates) would be very beneficial. This system is used by two applied maths communities: "finite elements computing" and - increasingly - "computer graphism".
It means converting the cartesian formulation (x,y,z) into a formulation with (a,b,c,d) (normalized by $a+b+c+d=1$ that can be given a physical interpretation as "weights").
What is their principle? Being given a point M inside (or outside) ABCD, its barycentric (or tetrahedral) coordinates are the (unique!) (a,b,c,d)=(a,b,c,1-a-b-c) such that M is the weighted average of $A$ with weight $a$, $B$ with weight $b$... This is written: M=aA+bB+cC+dD;
This could be done either by fixing a reference tetrahedron and computing all with respect to it, or define a system of barycentric coordinates with respect to each tetrahedron of your "tetrahedrization" ;
Barycentric coordinates behave well with respect to linearity and convexity.
3 references among many:
A nice presentation that can be found on the internet :
A lot of pointers can be found in : https://gist.github.com/paniq/8bdec20d00d08810f081
See also the article: "Barycentric coordinates for convex sets" by Warren, Schaefer, Hirani, Desbrun, in Advances in Computational Mathematics (2007) 27: 319-338