Tetrahedra of varying density and their applications

We propose concepts to utilize basic mathematical principles for computing the exact mass properties of objects with varying densities. For objects given as 3D triangle meshes, the method is analytically accurate and at the same time faster than any established approximation method. Our concept is based on tetrahedra as underlying primitives, which allows for the object’s actual mesh surface to be incorporated in the computation. The density within a tetrahedron is allowed to vary linearly, i.e., arbitrary density fields can be approximated by specifying the density at all vertices of a tetrahedral mesh. Involved integrals are formulated in closed form and can be evaluated by simple, easily parallelized, vector-matrix multiplications. The ability to compute exact masses and centroids for objects of varying density enables novel or more exact solutions to several interesting problems: besides the accurate analysis of objects under given density fields, this includes the synthesis of parameterized density functions for the make-it-stand challenge or manufacturing of objects with controlled rotational inertia. In addition, based on the tetrahedralization of Voronoi cells we introduce a precise method to solve L2|∞\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_{2|\infty }$$\end{document} Lloyd relaxations by exact integration of the Chebyshev norm. In the context of additive manufacturing research, objects of varying density are a prominent topic. However, current state-of-the-art algorithms are still based on voxelizations, which produce rather crude approximations of masses and mass centers of 3D objects. Many existing frameworks will benefit by replacing approximations with fast and exact calculations.


Introduction
When it comes to estimating masses or computing gravitational centers for objects of varying density almost everywhere discretized approximations are used. Often their accuracy is not even questioned. Surprisingly, it is neither a very complex problem nor bound to numerical tradeoffs between computation cost and accuracy. As used in finite element methods (FEM) [23,24], our techniques are based on tetrahedra as the underlying geometry primitive. Density fields are accurately represented by specifying the density at the four vertices. The analytically derived closedform integrals for mass properties are expressed with simple matrix-vector multiplications. In Sect. 2, we derive the concept, demonstrate its versatile applicability in Sect. 3 and conclude in Sect. 4 with numerical comparisons and a quantitative evaluation.

Contributions
In our paper, we briefly review the mathematical basis that allows for accurate solutions of mass and mass-center integrals of objects with varying density. It is simple as computations for tetrahedra with constant density are easily extended towards varying density. It is precise due to the use of a tet-mesh itself as input and therefore avoiding aliasing bias from discretization, as it is common in stateof-the-art applications using axis-aligned voxelizations. It is fast as computations can be implemented as matrix-vector multiplications, suitable for vectorized or GPU execution. It is versatile because the framework is not limited to density and can be extended to integrate arbitrary linear properties over a volume.
Our application framework, introduced in Sect. 3, can be summarized with the following main contributions: • Arbitrary Objects The solutions, known for tetrahedra, can be straightforward generalized to arbitrary polyhe-dra. We can accurately determine mass properties for any tetrahedralized input object in specified density fields. • Optimizing Density We can also invert the problem and optimize a parameterized density field for an object and given mass properties. • Approaching nonlinearity Arbitrary nonlinear functions can be approximated in a Taylor-like piecewise linear fashion, limited in accuracy only by the tetrahedralization's resolution. • Additive manufacturing Combining the former two contributions, we eventually propose to 3D-print objects of nonlinear density with optimized mass properties for balance or rotation-aligned inertia momentum. • Expressing energy functions As the method is not bound to only physical characteristics, we extend a Lloyd relaxation procedure based on L 2 Voronoi cells, using the concept to replace the cells energy integral with an accurate closed form solution for an L ∞ objective function.

Related work
Over the last few years, 3D printing not only attracted the doit-yourself hobbyist community but also gained popularity in various industrial applications. Nowadays, additive manufacturing processes go way beyond stacked layers of plastic and support a wide range of multiple or mixed materials, even including metals. Its widespread use, e.g., in the medical [31] or automotive [11] industry, keep this a relevant research topic. Hence the general interest of the computer graphics community for analyzing and processing 3D geometry, research in this field also spawned state-of-the-art algorithms aiming at 3D manufacturing. The procedure introduced by Prévost et al. [22] allows 3D models to be balanced in a specific position by shifting the object's center of gravity over a safe area on which the object eventually stands. Optimized weight distribution is achieved by carving out the object's interior and deformations of the hull, if necessary. Advancements of this technique optimize objects to have rotation-symmetric weight distributions and allow them to spin-like toy-tops [1]. Multistable balancing states are accomplished by using movable masses [21]. However, established techniques for mass property optimization are still based on approximations of the actual volume and mass distribution using quantized voxelizations. The approach by Musialski et al. [18] utilizes offset surfaces for shape and mass property optimization but also relies on binary material distribution. Even publications specialized on varying density for manufacturing [10,29] discretize their density field with marching cubes [15] or octrees [16] combined with dithering techniques. Known methods for computing exact mass properties of polyhedral bodies [17] are restricted to constant density. Like ours, they are based on integrals over the volume and surface of an object. A later revision [6] made the concept feasible for implementation. The varying density of polyhedral bodies, however, was first studied in the field of geophysics and concluded with the focus on gravitational fields [5,7] but not general mass properties. Our approach to computing accurate masses and mass centers under varying density relies on a tetrahedral decomposition of the input object. TetGen [26] is a tetrahedralization tool for polyhedral manifolds based on a Voronoi/Delaunay tessellation. Most recently, TetWild [9] introduced another fast and robust way to tetrahedralize any given 3D triangle soup, providing many adjustable parameters to the user. Our results for examples and applications are based on the outputs of both tools but often specifically on TetGen since it is able to preserve the original input surface. However, any other tet-meshing pipeline will generate equally suitable input as well.
As mentioned with TetGen, tetrahedralization is closely related to Voronoi and Delaunay graphs. In our Sect. 3.5, we introduce a novel approach on Lloyd relaxations (based on Voronoi tessellations) using the L ∞ norm. Ray et al. proposed to compute meshless Voronoi [25] and restricted power diagrams [2] on the GPU, however, both only in the common L 2 space. It definitely is a promising task to explore combinations of their diagrams and our take on L ∞ relaxations.

Concept
Our goal is to approach mass properties for polyhedral manifolds of varying densities with analytical tools. In order to compute the mass, the center of mass, or other related quantities in a field of varying density that is bounded by a triangle mesh we use closed-form solutions for a tetrahedralization, induced by the given mesh surface. These kinds of approaches are admittedly standard in FEM but so far rarely have been used in computer graphics. Therefore, this section briefly summarizes all important formulas and introduces the geometric concept that allows for computing these quantities exactly for linearly varying density fields inside a tetrahedron. Detailed derivations of the resulting equations are featured in "Appendix A".

Problem statement
Mesh data structures are the straightforward and, therefore, most common way to store and represent 3D manifolds. Under constant density, mass properties like the center of mass or an inertia momentum can be easily computed directly on a mesh using the divergence theorem. With varying density, however, a volume bounded by arbitrarily shaped polyhedra cannot be directly integrated. A common fallback solution is to approximate the object shape by decomposi-

Fig. 2
A general tetrahedron T with varying density, defined by its four vertices v i . In addition to their geometrical dimensions (x, y, z), each vertex is attributed with a fourth density dimension d tion into feasible volume entities, usually voxels. However, as illustrated with the 2D example in Fig. 1, the success and precision of the approximation will always be limited by the chosen resolution (for space and values), even hierarchical concepts can only reduce sampling artifacts but not fully avoid them. In contrast, our solution for the computation under varying density is based on the simple idea of an alternative volume representation, namely the tetrahedron. As illustrated on the right in Fig. 1 with a trivially triangulated 2D shape, every 3D shape with a polygonal surface can be decomposed using tetrahedra. With tetrahedra, the mesh's true shape can be used in all computations and, therefore, corresponding results are free of discretization and aliasing bias.

Geometry integration
Computing mass properties for the general tetrahedron T , specified in Fig. 2, is trivial for constant density: For example with d i = 1, the mass is equal to the tet-volume and the center of mass is equal to its centroid. However, as the density attributes at each vertex can be individually specified, expressing a linear density field inside the tetrahedron, the computation of mass and mass center changes. Instead, as in FEM [23,24], we utilize a simple basis case in a linear density field, for which the integration is solved analytically. A linear combination of four base cases (one per vertex) already gives the desired properties for a general tetrahedron.

Mass properties for arbitrary tetrahedra
Mass Due to linearity, the mass for a tetrahedron with four different density values at its vertices can be simply expressed as the mean of these values times the volume, formulated in Eq. 1. (1) Center of Mass As expressed in Eq. 2, the mass center computes as a weighted sum of vertex positions and their normalized density values. An extended derivation for the combination of the four base cases to this general form is included in Eq. 12.

Application
With the presented approach, one can now calculate the mass and further mass properties of tet-meshes with arbitrary density fields efficiently and exactly. The power of the approach will be demonstrated in three different application scenarios: For example, as a fast and accurate replacement for widespread voxel-based approximations of arbitrary objects' mass properties. Further, it can be used to optimize the density distribution inside an object to obtain a given center of mass or a stable rotation axis, potentially even with nonlinear density fields. By introducing a closed-form solution, our concept even allows us to formulate an objective function in a volumetric Lloyd relaxation process, which was so far not analytically feasible.

Mass properties of arbitrary objects
With the techniques, introduced in Sect. 2, to compute mass and center of mass for general tetrahedra, we can generalize this concept further and approach arbitrary polyhedral manifolds: Objects are partitioned into tetrahedra, mass properties are determined individually, and results eventually recombined. Any tetrahedral mesh is suitable as input for our method; if the model is not already available as tetrahedral mesh, it can easily be generated using freely available tools like TetGen [26] or TetWild [9]. "Appendix C" proves the concept to be invariant of the actual tetrahedralization.
For an object O and a given density-field, one can now compute the accurate mass M T i and center of mass C T i for all tetrahedra T i ∈ O using Eqs. 1 and 2, respectively. These calculations can be executed very efficiently, using fast matrix-vector multiplications. As formulated in Eq. 3, the object's overall mass is obtained by simple summation and the center of mass as mass-weighted dot-product. Further, one may extend the derivation, as described by Tonon [28], for the inertia tensor T i of a general tetrahedron with specified density values of the individual vertices. Following the rules for rigid bodies and the parallel axis theorem, one can derive further mass properties, e.g., the moment of inertia as formulated with the inertia tensor O in Eq. 4, where I 3 is the 3×3 identity matrix and ⊗ the outer product.

Optimizing density fields
Now, that an object's center of mass can be determined for a given density field, one can invert the problem and fit a density field to an object where the mass properties are given.
As an exemplary use case, we approached the make-it-stand challenge described by Prévost et al. [22] to balance objects in a given pose. The center of mass has to be within certain boundaries of a projected surface polygon on which the object is supposed to be balanced. However, a solution to this problem is limited by the following constraints: (i) Negative mass is reasonable only in theoretical fields of physics, so we limit our model to the realm of positive density for now.
(ii) Zero density is a special case that can be modeled with our concept, e.g., with d i = 0. (iii) The shape of an object together with constraints (i) and (ii) will put some limits on the achievable location of an object's center of mass, e.g., it simply cannot pass a certain point. Rather than optimizing the per-vertex density directly, let's first consider a simplified density field as illustrated in Fig. 4. We utilize two planes that separate volumes of constant minimum d min and maximum density d max respectively, sandwiching a slice of volume of width r with linearly To simplify many computation steps we fix the density field to be axis aligned, i.e., the planes are parallel to the yz-plane. As accommodation for this fixed orientation of the field, the optimization needs to rotate and to translate the object accordingly instead. To realize the arbitrary location of the bisection-planes in the density function, we have to prepare our input mesh by intersecting some of the tetrahedra, see "Appendix B" for details. The energy to be minimized is formulated as the Euclidean distance E = |C − C O | between a target point C and the object's current center of mass C O when embedded in the density field, obviously with respect to the object's rotated and translated state.
The optimization energy is smooth and in some sense probably differentiable but deriving gradients is left for future work. In our experiments, we used Powell's method [20] to minimize the objective function. Figure 5 compares our balanced objects to crosssections of Prévost et al. Their proposed method found a solution to make the three spheres stand, by carving out the voxelized interior and deforming the object. To move the mass center of the Spheres into the balance region, the top sphere is shrunk and the bottom sphere is enlarged. For the MrHumpty figure to stand upside down, the belly is enlarged and half the interior carved out to compensate for the offaxis legs. In our results, the objects remain untouched as they are only embedded in an optimized density field. As our output geometry incorporates the input, error measures like the Hausdorff distance are simply 0.

Results
Our first experiment meets the same conditions as Prévost et al. where the center of gravity only has to be on the central vertical axis of the bottom sphere so that the object is in balance. For our next experiment we chose the center of mass to also be located centered in the bottom sphere, but 10% of the sphere's radius below its horizontal equator-line. Due to this low center of gravity, the standing spheres would roll into this position on their own (Table 1).
Our optimization managed to define density fields for which the object's center of mass is exactly on the specified axis or target point respectively. The results have regions of constant minimum and maximum density with a tilted and shifted gradient between them. Due to the symmetry of the  Fig. 5 An unstable input of three spheres and a figure which is supposed to stand upside down. Prévost et al. [22] managed to balance the objects by deforming them, caving out the interior and shifting the center of mass over a defined safe-region. Our balanced version of the spheres can also be balanced on a small flattened face, the standing version will roll into this position on its own, due to its low center of gravity. Varying density is sufficient to balance the objects, deformation is not required. As reference for Table 2, the spheres are scaled to have a diameter of 1 and MrHumpty has a hand-to-hand width of 2  objects, the angle β is zero. To approach somewhat reasonable manufacturing limits, we set t = 1 (= d min ). The other found parameters are given in Table 2. Results of this comparison should be seen as a theoretical proof of concept, as this rather unconstrained optimization leads to quite high values for the gradient steepness s. Density differences of this multitude are ill-suited for current single-material manufacturing techniques. Additive multi-material techniques, on the other hand, could approximate smooth gradients like this, e.g., using dithering. The rockerarm is a prominent example for an asymmetric object with rotary mount. Due to imbalance, the native center of mass (red) is not located on the rotation axis. With optimized density, our center of mass (green) is located on the rotation axis and the principle inertia momentum axis is parallel to the rotation axis core with a Gaussian slope. Tet-mesh vertices become 3D sampling positions for the 3D density field, however, gradients within each tetrahedron are still linear. Nevertheless, this piecewise linear Taylor approximation of a nonlinear field is C 0 continuous everywhere (C ∞ within a tetrahedron). The accuracy of this representation is only limited by the resolution of the tetrahedral mesh, which can be specified in common tetrahedralization tools.

Optimizing nonlinear fields
Advanced applications, specifying more than a single center of mass, may require density fields, more sophisticated than linear gradients. An approach related to the makeit-stand concept proposed the challenge to make objects spinnable [1] by moving mass centers to a specified rotation axis. This is not only a desirable criterion for toy tops or yo-yos but is also of great value in any mechanical process involving rotating movements to reduce the wear and tear of involved components. Engineering such mechanical components often comes with tight constraints on available space and does not allow for arbitrary placement of counterweights. Figure 7 illustrates an example with the rockerarm object, which is to be mounted on rotary bearings. With constant density, the native center of gravity and inertia tensor are off-axis due to the obvious asymmetry of the object.
For this object, the optimized density field results in a center of mass located on the rotation axis along with a parallel principle inertia momentum axis. The density field is parameterized with the nonlinear density function d(v) (Eq. 5), where v is a tet-mesh vertex, p a 3D coordinate and k a scalar factor.

Results
Optimizing for a specific center of mass, as in Sect. 3.2, is not trivial but possible, dependent on given constraints. Additionally fitting a principle inertia momentum axis, however, can be challenging as the density distribution for a certain center of mass may be in conflict with the optimal density for the momentum axis. A field parameterization with more degrees of freedom than ours (Eq. 5) might be more suitable for optimization but unreasonable for practical results. Our results are shown in Fig. 7 with an optimized center of mass (green). The density parameters are: The principle inertia momentum axis was met with accuracy of < 1 • , the center of mass is actually precisely located on the specified target rotation axis.

3D printing of varying density
For some objects, we demonstrate both synthetic results as well as 3D printouts. One has to mention that 3D printing hardware for varying density is still in an early development state [8,13,19] and the range of available densities is limited. On recent Prusa FDM printers, however, it is possible to alter the extrusion rate while printing. The first step is to obtain the G-code for an object with a regular infill pattern. In order to approximate the optimized density field, we modify the Fig. 8 Modifying the extrusion rate allows for printing density gradients. The bar on top is a cross-section of an object 10×10×100 mm in size, scanned with a photocopier. Pictures of the rockerarms (62.9×33×17.6 mm) were taken in front of a lightsource to highlight the different density distributions. See video for live demo line thickness to vary along printed segments by accordingly adjusting the relative extrusion rate in the G-code slice by slice. Figure 8 shows a 3D printed example of a simple bar with increasing density. The varying amount of printed material alters the translucency of the object. The 3D printed rockeram of Fig. 7 with optimized density for an on-axis center of mass and a parallel principle inertia momentum axis leads to significantly smoother spinning as can be seen in the accompanying video.

Lloyd relaxation with the L ∞ norm
The last proposed application scenario makes use of the closed-form integral solution of measures on volumes of varying density to approach Lloyd relaxations under nonstandard norms. All calculations are done on the actual shape of the Voronoi cells avoiding any voxelization, which reduces artifacts and speeds up the computation. Lloyd's algorithm [14] is an iterative optimization procedure that is proven to converge to Centroidal-Voronoi-Tessellations (CVT) under the L 2 norm [4]. The iteration alternates two steps: I. Compute a Voronoi diagram for a given set of points. II. Reposition each point to the centroid of its Voronoi cell. This can also be formulated as an optimization task, minimizing the diagram's global energy function.
Since the native L 2 cells are all convex, the computation of new centroids is quite simple. However, in many meshing applications, L p or even L ∞ are more desirable [12], due to their more rectangular or cubical cell shapes. For L p norms ( p > 2), the Voronoi cells are no longer always convex and the diagram becomes very impractical to handle or even generate since there is (to the best of our knowledge) no software library that is able to compute L p or L ∞ Voronoi tessellations.
In meshing applications [12,27] the diagram itself is actually not relevant, but only the site positions are of interest [25]. We propose a way to compute Lloyd relaxed site positions with the L ∞ metric, also called the Chebyshev distance: First, cell geometry and topology are borrowed from an L 2 Fig. 9 2D visualization of the equivalent energy terms of Eq. 7 with an integral over the Chebyshev distances of all points in a cell or the sum over its tetrahedral fragments with density gradients   Left: Average movement of all cell centers in one iteration, given in percent of the optimal cell's diagonal. Right: Average L ∞ energy of all cells, given in percent of an optimal cell's L ∞ energy, which is why the result converges to 100% tessellated diagram, which comes with the convenience of convex-only cell shapes. Then we use our method and compute a cell's mass, reinterpreted as the energy which is to be minimized by a new cell center.

The Energy Term
For the goal to minimize the L ∞ energy within a cell, let us briefly recapitulate how the Chebyshev distance d ∞ is defined. As formulated in Eq. 6, the distance between two points p and q is the maximum of their absolute differences over all dimensions, in the 3D case x, y and z: Equation 7 formulates the energy E C of a cell C as the total Chebyshev distance of all points P ∈ C to the cell's centroid C C . However, there are infinitely many points P ∈ C, so the energy function can only be evaluated with a nontrivial integral over the cell volume.
As illustrated in Fig. 9 (in 2D), this integral becomes feasible as a finite sum of analytical solutions. To achieve this, a cell is split into six fragments F k (k ∈ [±x, ±y, ±z]), as illustrated in Fig. 11. This effectively separates all points P within the cell with respect to their maximum differencedimension (Chebyshev). Due to the separation into the six fragments, the d ∞ (C C , P) distance dimension conveniently coincides with the corresponding geometric dimension k, i.e., the distance linearly increases along one of the coordinate axes. As the hull of a Voronoi cell might be complex, the six fragments are tetrahedralized using a trivial triangulation of their hull faces and the cell center itself. The inner sum in Eq. 7 accumulates masses M k T of all tetrahedra T in a  fragment F k as defined in Eq. 1. The Chebyshev distance is simply encoded as the density dimension along the coordinate axes for our computation. The outer sum accumulates the density-(or Chebyshev distance-) weighted volumes of the six fragments, resulting in the cell's L ∞ energy. The cell center is finally repositioned to minimize the computed L ∞ energy using the L-BFGS-B algorithm [3,32] for bound constrained minimization.

Results
Although cell energies are only optimized on an individual basis, the relaxation iteration also leads to a global decrease of the diagram's energy, analogously to the L 2 case. With our reformulation of the objective function, the second part of the Lloyd relaxation (repositioning of cell centers) becomes feasible for the L ∞ norm. The initialization of each iteration is still based on a computable L 2 Voronoi tessellation, which turned out to be sufficient as the relaxation still converges. Considering the shape of an L 2 Voronoi cell while optimizing the centers for the L ∞ energy, this optimization is not a full L ∞ relaxation but a convenient alternative. If an L ∞ tessellation was available, the relaxation would probably converge even faster and would also allow for individually oriented cells. Nevertheless, considering the alternatives, e.g., labeling underlying high-resolution voxel grids, it is an improvement in terms of both accuracy and performance. The plots in Fig. 12 show convergence results of the cube (Fig. 10b) throughout 50 Lloyd relaxation steps. The Movement plot shows the average distance traveled by all sites (cell centroids) in the Voronoi diagram during each relaxation step. This distance is given in percent of an optimal cubical cell's diagonal. The average cell-energy, shown in the Energy plot, is computed for each cell as illustrated in Fig. 11 by separating a cell into six fragments and accumulating the density-weighted volume. It is expressed in percent of an optimal cubical cell's L ∞ energy. Therefore, the convergence towards 100% indicates that our cells approach the anticipated optimal cubical cell shape. Fig. 13 By increasing the depth of an octree (# levels on the x-axis), the approximations converge against our analytic results. The plot shows the time in seconds to build and traverse the tree, the number of nodes in the tree, the mass center error C , the mass error M and inertia tensor error . Dashed lines show timings for our computation based on tets and the time it took for TetGen [26] to tetrahedralize the input

Discussion
This section presents the results of the proposed application scenarios for our concept. Benefits over traditional methods in terms of performance and accuracy are quantitatively discussed with numerical results. After an outlook on potential extensions and future work, we conclude with a roundup of our main contributions.

Results
As mentioned in Sect. 1.2, voxelizations or octrees are currently the most common method to approximate mass properties for objects in fields of varying density. Table 3 documents comparisons of our exact tetrahedron-based method to octrees of different depths. Our results provide the ground truth reference, to which the octree approximations are compared. Timings for the octrees include the build-up phase and traversal to compute the results. The most demanding parts in the build-up are in/out-tests, to decide if a cell is to be split again. To be comparable, we included the time to create the tet-mesh inputs for our method from basic surface meshes. Delaunay tetrahedralizations are computed with TetGen [26], using the -Y option which preserves the source mesh, so that our method and the octree have the exact same input. Both, octree and our method, are implemented in Python using vectorized NumPy arrays where possible for optimized efficiency. Although our method is well suited to be implemented in parallelized GPU code, all timings are measured on a single CPU core. The measurements show that, not only compared to the very deep but also for the small octree of only 4 levels, our method is multitudes faster even including the input tetrahedralization.
Timings and performance aside, the probably most valuable takeaway is the analytical accuracy of the results. Our method establishes actual ground truth results for mass properties under varying densities. Figure 13 plots the mass-property-errors of octree approximations (Fig. 14) converging against our results, as we increase their depth and therefore accuracy. Featured errors of mass M , center of mass C , and the inertia tensor are specified in Eq. 8.
In theory, a voxel grid of infinite resolution or an octree of infinite depth would give correct and unbiased mass property results. We use this capacity to show that octree results converge against our analytic results by increasing their depth and accuracy. The error plots do not converge monotonically due to aliasing artifacts, for some lower levels the approximations are just more accurate by chance.

Conclusion
We propose novel application scenarios for object's mass properties under varying densities. Easy to use analytical solutions make approximations obsolete, which are still common in recent state-of-the-art applications [1,10,21,22,29]. Our concept is fast, lightweight, easy to implement, and suitable for vectorized or parallelized frameworks. We demonstrate possible use cases where our method can be utilized straightforward: Masses, mass centers, and inertia tensors of arbitrary manifolds in given density fields are computed accurately. We formulate an optimization to determine a parameterized density field for an object and specified mass properties like a center of gravity or inertia tensor. Our proposed modification of the Lloyd relaxation is a novel L 2|∞ hybrid that allows us to imitate real L ∞ relaxations, which is a leap forward compared to the existing approximative alternatives. While our approach may find direct application in established research topics as meshing, spatial tessellation and simulation [9,25], we also see great potential in the young scientific field of additive manufacturing and hope to inspire many further research [10,29,31].
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
The authors declare that they have no conflict of interest.

A In-depth derivations
This appendix completes the derivations of the closed-form equations for mass and center of mass calculations used in Sect. 2. Figure 3 illustrates an exemplary basis case with density d D = 1 and 0 at the three other vertices. Since the density is normalized, it is not surprising that the density integral over the volume in Eq. 9 results in a quarter of the tetrahedron's volume V T . Thus, when combined for a general tetrahedron embedded in an arbitrary density field, the mass M T computes as the tetrahedron volume times the mean of all four density values, as formulated in Eq. 10.
Equation 11 formulates the volume integration of the mass center C T (D) for the shown base case, using a scaled center vector − → w .
The center of mass C T in Equation 12 for the general tetrahedron computes as the combination of the individual base case mass centers, weighted and normalized by the individual density values at the four vertices, respectively.

B Tet-split
To model the density function described in Sect. 3.2, the object is bisected using split-planes. As the input is already tetrahedralized, there might be tetrahedra that are only partially in one or the other region. This potentially violates the constraints of Sect. 3.2, e.g., if a tetrahedron would cross the 0-density limit. Therefore, such tetrahedra with vertices on both sides of a split-plane have to be cut. Figure 15 illustrates the four possible cases, how a plane may intersect a tetrahedron. The splits result in geometry that is always further tetrahedralizable, hence triangular prisms, quad-based pyramids or trivial tetrahedra. Affected tetrahedra in the input are easily identified by checking if they match one of these cases.
After the cut, corresponding tetrahedra are simply replaced by the subdivided geometry.

C Proof of concept
This appendix aims to demonstrate the validity of our approach, especially proving mass properties of polyhedra to be invariant of the used tetrahedralization. Therefore, we assemble the simple scenario shown in Fig. 16: An axisaligned box is (w.l.o.g.) centered on the x-axis with the AE H D quad face parallel to the yz-plane and the B FGC quad at a distance h to the origin at x 0 . The density gradually increases over the x dimension dependent on the density function d(x) = s · x + t. Box-related equations are denoted with an overset box-symbol , the tetrahedra equivalents are marked with a triangle-symbol .

C.1 Setup
For the constructed box scenario, it is fairly easy to determine its mass and mass center via integration as formulated in Eqs. 13 and 14, respectively. The cross-section-area of the box is formally expressed as a function B (x) over the integration domain, which is constant and, for simplicity, assumed to be 1.

C.2 Mass
To demonstrate equality of our approach to the box-solution, we first gather all the tetrahedra-related components. Equation 15 lists the independent masses of the tetrahedra, using the assumption of the box's cross-section to have an area of 1. The volume of the individual tetrahedra in this configuration compute as V rgcm = 1 6 and V y = 1 3 . They scale linearly if Masses of the individual tetrahedra (Eq. 15) summed up (Eq. 16) prove equality to the integrated box-mass (Eq. 13).

Center of mass
In Eq. 14 the center of mass integral results as a scaled vector − → w , parallel to the x-axis and scaled dependent on the density function. However, the center of mass for a general tetrahedron is formulated in Eq. 2 solely based on its vertices. To eventually express the equality of the box-and tetrahedralized solutions, we first reformulate the center of mass for tetrahedra in a similar − → w vector-dependent way. Therefore, we establish the three basic integrable tetrahedra cases shown in Fig. 17.