Computation of Implicit Representation of Volumetric Shells with Predefined Thickness

Abstract

We propose and validate a method to find an implicit representation of a surface placed at a distance h from another implicit surface. With two such surfaces on either side of the original surface, a volumetric shell of predefined thickness can be obtained. The usability of the proposed method is demonstrated through providing solid models of triply periodic minimal surface (TPMS) geometries with a predefined constant and variable thickness. The method has an adjustable order of convergence. If applied to surfaces with spatially varying thickness, the convergence order is limited to second order. This accuracy is still substantially higher than the accuracy of any contemporary 3D printer that could benefit from the function as an infill volume for shells with predefined thicknesses.

1. Introduction

Implicit surfaces are of considerable interest in engineering as they provide an efficient means for the modeling of complex shapes. However, if transferred to the physical world by any means of production, any surface must have a finite thickness and hence, becomes a volume. This can be done by tools for computer aided design by finding surface points and extruding them in the normal direction. For example, Payne and Toga [1] proposed an algorithm to compute a distance field for a discretized surface such that the offset surface could be extracted as an isosurface of the distance field. Venkatesh applied the marching cube algorithm to generate volume models from triply periodic minimal surfaces [2]. Held et al. [3] used generalized Voronoi diagrams to compute offset functions of variable thickness. Pham presented an overview on offsetting methods [4], which was later extended by Maekawa [5]. However, these reviews do not address the analytic offsetting of implicit surfaces.

Fayolle et al. [6] proposed a method for offsetting implicit functions by computing level sets of the distance function. The distance function is obtained by normalizing the original implicit surface by Rvachev’s normalization [7]. Fayolle’s method is based on numerical procedures that require sampling of the initial implicit surface. They apply the Euler forward method to propagate the offset into the direction of the normal to the surface.

Sethian [8] provides a good overview over fast marching methods used for offsetting implicit surfaces. He also discusses how the accuracy of the the fast marching method can be improved by discretizing the required normal derivatives with finite differences of higher order. However, this by itself does not guarantee an improvement of the order of convergence of the method.

In the mentioned methods, the resulting new surface is either a discretized surface with an explicit surface description or a discretized form of the implicit surface. In both cases, many of the appealing properties of the original implicit surface are lost.

Explicit surface representations usually require much more data than implicit representations, especially if the resolution requirements are high. Applications for implicit surfaces are in 3D printing [9,10,11] and simulation software. In both domains, tessellated surfaces are currently the state of the art. However, simulation software can suffer from artifacts when curved surfaces are approximated by triangles. For example, the accurate simulation of flow around a sphere at a high Reynolds number can be qualitatively wrong if the sphere is approximated by triangles due to spurious flow separation at the edges of triangles, while accurate results can be obtained with the sphere defined as geometric primitive [12], which also requires less memory.

The purpose of the current paper is to propose a method to obtain a new implicit surface at a predefined distance to a given implicit surface.

The reminder of the paper is organized as follows. In Section 2, we introduce a low-order approximation of the offset function to the implicit surface. In Section 3, it is shown how the approximation of the low-order scheme can be improved by Richardson extrapolation in order to obtain a second-order convergent method. In Section 4, it is shown how this methodology is extended to obtain offset functions at any order of convergence. Section 5 demonstrates empirically that the predicted convergence orders are in fact obtained by computing the error norm for offset functions for three different triply periodic minimal surfaces. In Section 6, we discuss the case of a variable thickness of the shell and demonstrate second-order of convergence. Conclusions are given in Section 7.

2. First-Order Offset Function

Let $f\left(\overrightarrow{x}\right)=0$ be an implicit surface $\Gamma$. Our aim is to find an implicit surface ${\Gamma }_h$ as the solution of the offset function $\varphi (\overrightarrow{x},h)=0$ such that the minimal Euclidean distance between any point on $\Gamma$ to ${\Gamma }_h$ is h. Further, we require that the distance between ${\Gamma }_h$ and ${\Gamma }_{-h}$ is $2h$. An exact solution to this problem is difficult to find in general and might not even exist, but we can start from a low-order approximation by computing the Taylor expansion of $f\left(\overrightarrow{x}\right)$:

$$\varphi (\overrightarrow{x},h)=f\left(\overrightarrow{x}\right)-h{\partial }_{n}f\left(\overrightarrow{x}\right)+\mathcal{O}\left(h^2\right),$$

(1)

where the index n denotes the normal to $\Gamma$, which is approximated to lowest order by the direction of the gradient of $f\left(\overrightarrow{x}\right)$ such that Equation (1) becomes

$$\varphi (\overrightarrow{x},h)=f\left(\overrightarrow{x}\right)-h\sqrt{(\nabla f\left(\overrightarrow{x}\right))·(\nabla f\left(\overrightarrow{x}\right))}+\mathcal{O}\left(h^2\right).$$

(2)

Omitting the error term and rewriting $(\nabla f\left(\overrightarrow{x}\right))·(\nabla f\left(\overrightarrow{x}\right))={(\nabla f\left(\overrightarrow{x}\right))}^2$ gives us the first-order offset function ${\varphi }^{\left(1\right)}(\overrightarrow{x},h)$, which can be interpreted as an eikonal equation:

$${\varphi }^{\left(1\right)}(\overrightarrow{x},h)=f\left(\overrightarrow{x}\right)-h\sqrt{{(\nabla f\left(\overrightarrow{x}\right))}^2}.$$

(3)

From ${\varphi }^{\left(1\right)}(\overrightarrow{x},h)=0$, a first-order accurate approximation of the offset surface denoted by ${\Gamma }_h^{\left(1\right)}$ can be computed.

We note here that computing an offset function via the propagation of a level set is far from new [13]. In the following, we demonstrate that this simple method can be improved to any specified order of accuracy.

3. Second-Order Offset Function

A first-order accurate solution to the offset problem is only acceptable if h is very small. In what follows, we present a systematic approach to rise the order of accuracy starting from ${\varphi }^{\left(1\right)}(\overrightarrow{x},h)=0$.

Let us first note that since ${\varphi }^{\left(1\right)}(\overrightarrow{x},h)$ is an implicit surface itself, we are able to apply the same methodology again in order to obtain an offset function of ${\varphi }^{\left(1\right)}(\overrightarrow{x},h)$. In the following recursive computation of the offset function ${\varphi }_n^{\left(m\right)}(\overrightarrow{x},h)$, we will indicate the level of the recursion by subscript n and the order of approximation by superscript $\left(m\right)$. With this, it is, for example, possible to cover the distance h in two steps of equal size. To achieve this, we recursively apply Equation (3) in itself to obtain

$$\begin{array}{ccc}\hfill {\varphi }_2^{\left(1\right)}(\overrightarrow{x},h)& =& {\varphi }_1^{\left(1\right)}(\overrightarrow{x},h/2)-\frac{h}{2}\sqrt{{(\nabla {\varphi }_1^{\left(1\right)}(\overrightarrow{x},h/2))}^2}.\hfill \end{array}$$

(4)

This second iteration ${\varphi }_2^{\left(1\right)}(\overrightarrow{x},h)$ has the same order of accuracy as ${\varphi }^{\left(1\right)}(\overrightarrow{x},h)$, as can be seen by substituting ${\varphi }_1^{\left(1\right)}(\overrightarrow{x},h/2)=\varphi (\overrightarrow{x},h/h)+\mathcal{O}\left(h^2\right)$ and applying the triangle inequality while realizing that h is not a function of $\overrightarrow{x}$, i.e., $\sqrt{{(\nabla (\varphi (\overrightarrow{x},h/h)+\mathcal{O}\left(h^2\right)))}^2}\le \sqrt{(\nabla {\left(\varphi (\overrightarrow{x},h/h)\right)}^2}+\mathcal{O}\left(h^2\right)$.

While ${\varphi }_2^{\left(1\right)}(\overrightarrow{x},h)$ has the same order of accuracy as ${\varphi }^{\left(1\right)}(\overrightarrow{x},h)$, it is expected to be more accurate as it uses half the step size. Adding more steps would raise the accuracy further but the methodology would quickly become unpractical as the algebraic complexity of the implicit surface function grows with each step. Instead, it appears more promising to combine the obtained solutions through Richardson extrapolation [14]. Richardson extrapolation is a convergence acceleration method for sequences of approximations obtained with different smallness parameters (e.g., step sizes). It is applicable if the order of convergence with respect to the smallness parameter is either known or can be obtained directly from the series of solutions. For this, it must be possible to express the error of the method as a Taylor series in the smallness parameter. It is also necessary that the error does not depend on terms not related to the smallness parameter. Therefore, only results obtained by the same deterministic method can be combined in a Richardson extrapolation.

We recall the general equation for Richardson extrapolation given a solution $g^{\left(q\right)}\left(h\right)$ obtained with smallness parameter h and a solution $g^{\left(q\right)}\left(\beta h\right)$ obtained with a smallness parameter $\beta h$, both with the same method of order q:

$$g^{(q+1)}=\frac{g^{\left(q\right)}\left(h\right)-{\beta }^{-q}{g}^{\left(q\right)}\left(\beta h\right)}{1-{\beta }^{-q}}.$$

(5)

In the current case, the refinement factor is $\beta =1/2$ and the convergence order is $q=1$. Note that both methods have to be evaluated at the same final location, which in the current case, is after a distance h. However, the solution $g^{\left(q\right)}\left(h\right)={\varphi }_1^{\left(1\right)}(\overrightarrow{x},h)$ is obtained using a single step while the second solution $g^{\left(q\right)}(h/2)={\varphi }_2^{\left(1\right)}(\overrightarrow{x},h)$ uses two steps of size $h/2$. After some algebra, this yields

$${\varphi }^{\left(2\right)}(\overrightarrow{x},h)=2{\varphi }_2^{\left(1\right)}(\overrightarrow{x},h)-{\varphi }_1^{\left(1\right)}(\overrightarrow{x},h)=f\left(\overrightarrow{x}\right)-h\sqrt{{(\nabla {\varphi }_1^{\left(1\right)}(\overrightarrow{x},h/2))}^2}.$$

(6)

A second-order approximation of the surface ${\Gamma }_h^{\left(2\right)}$ is obtained from ${\varphi }^{\left(2\right)}(\overrightarrow{x},h)=0$.

4. Arbitrary-Order Offset Function

There are various systematic ways to rise the convergence order further. We propose here using an extended Richardson extrapolation strategy. For this, we briefly return to the theory behind the Richardson extrapolation and explain how the same method can be used to obtain results at arbitrary order.

The general assumption behind the Richardson extrapolation is that the approximating function ${\varphi }_n^{\left(1\right)}(\overrightarrow{x},h)$ can be written as

$${\varphi }_n^{\left(1\right)}(\overrightarrow{x},h)=\varphi (\overrightarrow{x},h)+\sum _{i=1}^{\infty }{\mathcal{E}}_i{\left(\frac{h}{n}\right)}^i.$$

(7)

Here, n is the number of equidistant steps or iterations to reach a distance h such that the step size is $h/n$. The constants ${\mathcal{E}}_i$ are independent of the step size. Richardson extrapolation ignores all but the leading-order error to get a better approximation of $\varphi (\overrightarrow{x},h)$ from the weighted combination:

$${\varphi }^{\left(2\right)}(\overrightarrow{x},h)=w_1{\varphi }_{n1}^{\left(1\right)}(\overrightarrow{x},h)+w_2{\varphi }_{n2}^{\left(1\right)}(\overrightarrow{x},h).$$

(8)

The weights have to be recovered by inserting Equation (7) into Equation (8) and equating coefficients with regards to $\varphi (\overrightarrow{x},h)$ and ${\mathcal{E}}_1$, where $\varphi (\overrightarrow{x},h)$ is to be recovered while the error term is to be eliminated, i.e.,

$$\begin{array}{ccc}\hfill \varphi (\overrightarrow{x},h):& 1=& w_1+w_2\hfill \\ \hfill {\mathcal{E}}_1:& 0=& w_1\frac{1}{n_1}+w_2\frac{1}{n_2}.\hfill \end{array}$$

(9)

Solving this set of equations gives rise to the usual Richardson extrapolation formula.

While Richardson extrapolation usually only considers the elimination of the leading order error term, the formula can be extended to cover several error terms at the same time. The conditions for the extended Richardson extrapolation can be written simply by

$$\begin{array}{ccc}\hfill \varphi (\overrightarrow{x},h):& 1=& \sum _{j=1}^m{w}_j\hfill \\ \hfill {\mathcal{E}}_i:& 0=& \sum _{j=1}^m\frac{w_j}{{\left(n_j\right)}^i},\hfill \end{array}$$

(10)

where $i\le m$. By solving Equation (10), the weights for an extended Richardson extrapolation at arbitrary order can be determined.

Table 1. Table of the weights for extended Richardson extrapolation up to order four.

Order	${\mathit{\varphi }}_1^{\left(1\right)}(\overrightarrow{\mathit{x}},\mathit{h})$	${\mathit{\varphi }}_2^{\left(1\right)}(\overrightarrow{\mathit{x}},\mathit{h})$	${\mathit{\varphi }}_3^{\left(1\right)}(\overrightarrow{\mathit{x}},\mathit{h})$	${\mathit{\varphi }}_4^{\left(1\right)}(\overrightarrow{\mathit{x}},\mathit{h})$
1	1	-	-	-
2	$-1$	2	-	-
3	$\frac{1}{2}$	$-4$	$\frac{9}{2}$	-
4	$-\frac{1}{6}$	4	$-\frac{27}{2}$	$\frac{32}{3}$

lists the weight coefficients for extended Richardson extrapolations up to order four.

As higher derivatives are successively eliminated from the Taylor expansion through the extended Richardson extrapolation, our method can also be interpreted as a way to obtain successively accurate approximations of $\sqrt{{(\nabla f\left(\overrightarrow{x}\right))}^2}$:

$$\sqrt{{(\nabla f\left(\overrightarrow{x}\right))}^2}=\frac{f\left(\overrightarrow{x}\right)-{\varphi }^{\left(n\right)}(\overrightarrow{x},h)}{h}+\mathcal{O}\left(h^n\right).$$

(11)

In this way, our method can be seen as a semianalytic way to construct successively more accurate approximations to the eikonal equation.

4.1. Implicit Infill Volumes

An implicit function for the infill volume can be derived from the implicit offset functions in positive and negative directions indicated by $h<0$ and $h>0$. We define the implicit volume $\Phi (\overrightarrow{x},h)$ as the shell of the zero contour of $f\left(\overrightarrow{x}\right)$ with thickness h by

$$\Phi (\overrightarrow{x},h)=\varphi (\overrightarrow{x},h/2)\varphi (\overrightarrow{x},-h/2).$$

(12)

The point $\overrightarrow{x}$ is said to be inside the shell if $\Phi (\overrightarrow{x},h)<0$ and outside the shell if $\Phi (\overrightarrow{x},h)>0$.

4.2. Limitations

The problem to supply an arbitrary surface with a finite thickness is not well posed in general. In particular, it is difficult to unambiguously define what is meant by thickness of a curved surface. For example, if the local radius of curvature of a surface is small compared to the desired thickness of the shell, the resulting body could be self-intersecting. In the proposed method, the thickness is defined through Equation (1) as the path length of a curve following the normal of the implicit function $f\left(\overrightarrow{x}\right)$. For shells with local radii much larger than the thickness h, this path is approximately straight, but for increasing thicknesses, the path itself can be curved and the distance of the end points of a curved path is always shorter than the length of the path. Hence, the method proposed here is only valid for shells with a thickness much smaller than their local radii of curvature. It is important to note that this limitation is intrinsic to the definition of thickness in the method and this limitation is not affected by the order of the Richardson extrapolation chosen. Further, we note that the method does not check for self-intersection.

5. Examples

In this section, we will investigate the performance of the proposed methodology considering some examples of varying complexity.

5.1. Circle

A circle can be defined in a Cartesian coordinate system by the implicit function

$$f(x,y)=x^2+y^2-R^2,$$

(13)

with R being the radius. For the zero contour of $f(x,y)$, it is straightforward to compute the two contours with offset h by replacing the radius with $(R-h)$ and $(R+h)$, respectively.

In order to validate the proposed method, we compute the first- and second-order offset functions from Equation (13):

$$\begin{array}{ccc}\hfill {\varphi }^{\left(1\right)}(x,y,h)& =& x^2+y^2-R^2-2h\sqrt{x^2+y^2}\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill {\varphi }^{\left(2\right)}(x,y,h)& =& x^2+y^2-R^2-2h\sqrt{x^2+y^2+\frac{h^2}{4}-h\sqrt{x^2+y^2}}.\hfill \end{array}$$

(15)

Rewriting Equations (14) and (15) in terms of polar coordinates gives

$$\begin{array}{ccc}\hfill {\varphi }^{\left(1\right)}(\rho ,\theta ,h)& =& {\rho }^2-R^2-2h\rho \hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill {\varphi }^{\left(2\right)}(\rho ,\theta ,h)& =& {\rho }^2-R^2-2h\sqrt{{\rho }^2+\frac{h^2}{4}-h\rho }.\hfill \end{array}$$

(17)

At the specified location of the new surface, the offset functions are

$$\begin{array}{ccc}\hfill {\varphi }^{\left(1\right)}(R+h,\theta ,h)& =& -h^2\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill {\varphi }^{\left(2\right)}(R+h,\theta ,h)& =& 0,\hfill \end{array}$$

(19)

which means that the second-order offset function provides the exact solution while the error of the first-order offset function is within the predicted range.

5.2. Triply Periodic Minimal Surfaces (TPMSs)

Periodic minimal surfaces have many potential applications in engineering. There are several ways to define a minimal surface [15]. It can be defined as the surface that locally has the smallest area. Another definition is that a minimal surface has a zero mean curvature. In the following examples, we assume that the TPMS is given by the zero contour of the function $\mathfrak{F}\left(\overrightarrow{x}\right)$. Without loss of generality, we assume the TPMS to be periodic in a cubic box with side length L.

The offset functions up to fourth order of convergence for a TPMS are obtained following the proposed procedure.

We compute four different first-order approximations of the offset function according to

$$\begin{array}{ccc}\hfill {\varphi }_1^{\left(1\right)}(\overrightarrow{x},h)& =& \mathfrak{F}\left(\overrightarrow{x}\right)-h|\nabla \mathfrak{F}\left(\overrightarrow{x}\right)|\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill {\varphi }_2^{\left(1\right)}(\overrightarrow{x},h/2)& =& {\varphi }^{\left(1\right)}(\overrightarrow{x},h/2)-\frac{h}{2}|\nabla {\varphi }^{\left(1\right)}(\overrightarrow{x},h/2)|\hfill \end{array}$$

(21)

$$\begin{array}{ccc}\hfill {\varphi }_3^{\left(1\right)}(\overrightarrow{x},h/3)& =& {\varphi }_2^{\left(1\right)}(\overrightarrow{x},h/3)-\frac{h}{3}|\nabla {\varphi }_2^{\left(1\right)}(\overrightarrow{x},h/3)|\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill {\varphi }_4^{\left(1\right)}(\overrightarrow{x},h/4)& =& {\varphi }_3^{\left(1\right)}(\overrightarrow{x},h/4)-\frac{h}{4}|\nabla {\varphi }_3^{\left(1\right)}(\overrightarrow{x},h/4)|.\hfill \end{array}$$

(23)

From this, we obtain different orders of approximation using the extended Richardson extrapolation:

$$\begin{array}{ccc}\hfill {\varphi }^{\left(1\right)}(\overrightarrow{x},h)& =& \mathfrak{F}\left(\overrightarrow{x}\right)-h|\nabla \mathfrak{F}\left(\overrightarrow{x}\right)|\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill {\varphi }^{\left(2\right)}(\overrightarrow{x},h)& =& 2{\varphi }_2^{\left(1\right)}(\overrightarrow{x},h/2)-{\varphi }_1^{\left(1\right)}(\overrightarrow{x},h)\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill {\varphi }^{\left(3\right)}(\overrightarrow{x},h)& =& \frac{9}{2}{\varphi }_3^{\left(1\right)}(\overrightarrow{x},h/3)-4{\varphi }_2^{\left(1\right)}(\overrightarrow{x},h/2)+\frac{1}{2}{\varphi }^{\left(1\right)}(\overrightarrow{x},h)\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill {\varphi }^{\left(4\right)}(\overrightarrow{x},h)& =& \frac{32}{3}{\varphi }_4^{\left(1\right)}(\overrightarrow{x},h/4)-\frac{27}{2}{\varphi }_3^{\left(1\right)}(\overrightarrow{x},h/3)+4{\varphi }_2^{\left(1\right)}(\overrightarrow{x},h/2)-\frac{1}{6}{\varphi }_1^{\left(1\right)}(\overrightarrow{x},h).\hfill \end{array}$$

(27)

The derivation of ${\varphi }^{\left(n\right)}(\overrightarrow{x},h)$ is easily automated in a computer algebra system such like Mathematica, Maple, or SymPy.

Mathematica source code for the generation of offset functions at different orders of convergence is presented in the Supplementary Material published alongside this paper.

5.2.1. Error Measurement Procedure

To measure the accuracy of the proposed method, we employ a sampling technique. We want to measure the distance between the surface of the TPMS and the zero contour of its offset function. The computer algebra software Mathematica is used to find a large sample of points N on one of the surfaces. The sampled points are uniformly distributed in x and y directions. The z coordinate is determined by solving for z given x and y on the surfaces. The procedure for collecting the points is described in Algorithm 1.

Algorithm 1: Points sampling from an implicit surface $F(x,y,z)=0$

input: implicit surface function $F(x,y,z)$ output: list of 3-tuple $SampledPoints$ //$\phantom{(}$ create two list of coordinates in x and y directions$\phantom{(}$ $A:=\{a\mid a\in [0,0.05,0.1,\dots ,0.95,1\left]\right\}$ $B:=\{b\mid b\in [0,0.05,0.1,\dots ,0.95,1\left]\right\}$ //$\phantom{(}$ substitute each a and b into F and solve for the c$\phantom{(}$ $SampledPoints:=\left\{\right(a,b,\mathrm{Solve}\left(F\right(a,b,c)=0,c))\mid a\in a\mathrm{and}b\in B\}$

From each point, a ray is shot in the normal direction to the surface. Then, the intersection distance $h_m^{\left(n\right)}$ with the other surface is measured and compared to the predefined distance h. The error is measured in the following norm:

$$\begin{array}{c}\hfill Error=\left|\right|h^{\left(n\right)}-h{\left|\right|}_2=\frac{1}{N}\sum _{i=1}^N\sqrt{{\left(\frac{(h_m^{\left(n\right)}\left(i\right)-h}{h}\right)}^2},\end{array}$$

(28)

where N is the number of points used and n is the order of the approximation used to represent the offset surface. The procedure for calculating the error is described in Algorithm 2. The error measurement procedure, including sampling of the points from the surface and finding the intersection distance to the next surface, is done on both surfaces, i.e., the original zero contour of the TPMS and the offset function.

Algorithm 2: Calculation of the $L^2$ error in distance h from points in the surface $F_1\left(\overrightarrow{x}\right)=0$ to a distanced surface $F_2\left(\overrightarrow{x}\right)=0$

5.2.2. Gyroid

A Gyroid is a TPMS that is described by the implicit equation

$$\begin{array}{ccc}\hfill {\mathfrak{F}}_G\left(\overrightarrow{x}\right)& =& sin\left(\frac{2\pi x}{L}\right)cos\left(\frac{2\pi y}{L}\right)+sin\left(\frac{2\pi y}{L}\right)cos\left(\frac{2\pi z}{L}\right)\hfill \\ & & +sin\left(\frac{2\pi z}{L}\right)cos\left(\frac{2\pi x}{L}\right)=0,\hfill \end{array}$$

(29)

where L is the edge length of a cube circumscribing a unit cell of the Gyroid shown in Figure 1 together with the other surfaces investigated in this paper. With the emergence of 3D printing, Gyroids are of great interest since Gyroid-based patterns can be used to fill space while providing sufficient strength and low weight to the 3D-printed parts. Figure 2 shows the bodies generated with the third-order offset function of the Gyroid for various thicknesses.

We employ the error measurement from Section 5.2.1 to validate the proposed method. Figure 3 shows the convergence of the different offset functions for varying thickness. It is observed that the predicted order of convergence is recovered. The first-order offset function appears to be converging with second-order for large enough thicknesses but approaches the predicted asymptote for thin shells.

5.2.3. Schwarz P

The next example is the Schwarz P TPMS, defined as

$${\mathfrak{F}}_{Sch}\left(\overrightarrow{x}\right)=cos\left(\frac{2\pi x}{L}\right)+cos\left(\frac{2\pi y}{L}\right)+cos\left(\frac{2\pi z}{L}\right)=0.$$

(30)

The resulting bodies obtained from the third-order offset function are shown in Figure 4. Again, we measure the convergence of the proposed method following Section 5.2.1 and display them in Figure 5. In this case, the asymptotic behavior of the four tested variants follows the theoretical prediction very closely.

5.3. Neovius

The previous examples had very homogeneous curvatures. For the next example, we chose a surface for which the curvature changes. The Neovius surface is defined as

$$\begin{array}{ccc}\hfill {\mathfrak{F}}_{Neo}\left(\overrightarrow{x}\right)& =& 3\left[cos\left(\frac{2\pi x}{L}\right)+cos\left(\frac{2\pi y}{L}\right)+cos\left(\frac{2\pi z}{L}\right)\right]\hfill \\ & & +4\left[cos\left(\frac{2\pi x}{L}\right)cos\left(\frac{2\pi y}{L}\right)cos\left(\frac{2\pi z}{L}\right)\right]=0.\hfill \end{array}$$

(31)

The bodies generated from the third-order offset function are depicted in Figure 6. It is seen that the intersection of the Neovius surface with the axes plains has an elliptical shape. Figure 6c shows how the hollow connections between unit cells are filled up. It is obviously difficult to guarantee a constant thickness of the shell in these connections once the thickness approaches the width of the openings. This is evident from the convergence analysis, executed in the same way as for the other two bodies and displayed in Figure 7. For small enough thicknesses, the theoretically predicted convergence is well captured. However, for thicknesses of $0.1$ L, the error deteriorates as the resulting shape can no longer be described as a shell of the Neovius surface. One observation is that this deterioration is less severe when the thickness is measured from the offset function to the Neovius surface than when it is measured the other way around. The reason for this is most likely seen in the fact that the proportion of the offset function contour with inaccurate distance towards the Neovius surfaces shrinks when the connections are closed as these closed areas are no longer part of the offset function.

6. Shells with Variable Thickness

It can sometimes be of interest to allow for a variable thickness of the shell [16,17]. Applying the method outlined in this paper to shells with variable thickness would imply substituting h by the function $h\left(\overrightarrow{x}\right)$ in all stages of ${\varphi }^{\left(n\right)}(\overrightarrow{x},h\left(\overrightarrow{x}\right))$ before taking the derivative. The implicit dependence of h on $\overrightarrow{x}$ would then automatically be taken into account. However, this approach would not produce the desired result, which can be explained as follows: The proposed method approximates the thickness of the shell by moving a predefined distance along the normal of $\varphi (\overrightarrow{x},h)=0$ for increasing h. This procedure is only accurate under the assumption that the normal of $\varphi (\overrightarrow{x},h)=0$ does not deviate substantially from the normal of $f\left(\overrightarrow{x}\right)=0$, such that the procedure follows a straight line. However, if h is not constant, the normal becomes, through the chain rule of differentiation:

$$\overrightarrow{n}=\frac{\nabla \varphi (\overrightarrow{x},h\left(\overrightarrow{x}\right))}{|\nabla \varphi (\overrightarrow{x},h\left(\overrightarrow{x}\right)\left)\right|}=\frac{\nabla \varphi (\overrightarrow{x},h)+\frac{\partial \varphi (\overrightarrow{x},h)}{\partial h}\nabla h\left(\overrightarrow{x}\right)}{|\nabla \varphi (\overrightarrow{x},h)+\frac{\partial \varphi (\overrightarrow{x},h)}{\partial h}\nabla h\left(\overrightarrow{x}\right)|}.$$

(32)

The direction of the normal will therefore be governed by how h changes with $\overrightarrow{x}$, which is another way of saying that if h is not constant, $\varphi (\overrightarrow{x},h\left(\overrightarrow{x}\right))=0$ is naturally not parallel to $f\left(\overrightarrow{x}\right)=0$.

Another and more straightforward approach to obtain variable thickness is to assume h to be constant in the derivation of $\varphi (\overrightarrow{x},h)$ and only replace it when evaluating $\varphi (\overrightarrow{x},h)=0$. This also has limited accuracy because if h is a function of $\overrightarrow{x}$, it changes by a value $\sim h{\partial }_{n}h$ over the thickness h. For smooth variations of the thickness, it is reasonable to assume that ${\partial }_{n}h=\mathcal{O}\left(h\right)$ such that $h{\partial }_{n}h=\mathcal{O}\left(h^2\right)$. Hence, the accuracy of the proposed method for variable thickness is bound to second order. This is still sufficiently accurate for many applications.

Figure 8 shows the convergence for shells with a linearly increasing thickness along the x axis for the different TPMS. It is observed that the first-order method remains to be first-order accurate, albeit it performs better than expected for the Gyroid at larger thicknesses. The third-order method typically shows better results than the second-order method but it too converges only with second order. Since the intrinsic error in the thickness is already of order $\mathcal{O}\left(h^2\right)$, further increasing the order of ${\varphi }^{\left(n\right)}(\overrightarrow{x},h)$ would not improve the asymptotic nature of the result. For a shell thickness of $0.01$ L, the second-order method already obtains errors smaller than $0.03\%$ of the target thickness in all three cases. Figure 9 shows the resulting bodies with variable thickness.

7. Conclusions

In this paper, we proposed a systematic methodology for deriving analytic implicit body shells of predefined thickness from implicit surfaces. The convergence order and thereby the accuracy of the thickness of the shell can be freely adjusted trough the extended Richardson extrapolation. The motivation for our method comes from engineering applications of additive manufacturing, where such infill volumes are required for the actual 3D-printing process on the one hand and for numerical simulations of such structures on the other hand. For constant thicknesses, our method reaches a precision of smaller than ${10}^{-6}$ for wall thicknesses measuring $1\%$ of the wave length of the periodic structure. This precision by far exceeds the accuracy of contemporary 3D-printing technology. By opting for offset functions of lower-order computational time can be reduced. The memory requirements for storing the offset function of a TPMS is negligible compared to a discretized surface mesh. Being analytic, the offset function does not lose any quality when scaled. At the time of this writing, it is not straightforward to use the offset function directly as input for commercial 3D printers as they typically expect a tessellated mesh representation of the surface, which can, of course, easily be obtained from discretizing the zero contour of the offset function. Our method can be used to identify, for any point in the domain, whether it lies in the shell or outside the shell. As the shell can be expressed as a solid body in this way, the question of water tightness of the surface never arises. However, the slicing of the object for G-code generation and 3D printing was not addressed in this theoretical work. In the presented semianalytic form, our method can only be applied to surfaces provided as differentiable implicit functions. Such surfaces are always closed as they provide an inner and an outer side. In 3D printing, it is often of interest to include open surfaces [18]. In our approach, this is possible by cutting the resultant shell with an additional solid in a postprocessing step.

Besides being useful for the definition of infill volumes for the printing process itself, the proposed method has other applications, for example, in the analysis and simulation of 3D-printed objects. In general, methods of adaptive manufacturing come with a finite smallest line width. In many applications, the line width of the printer is designed to coincide with the thickness of the shell structure under construction. This is, for example, common for 3D printing of slender concrete structures in applications of civil engineering and architecture [19,20]. But the same is also true for wire arc additive manufacturing [21] and the various types of powder bed methods like selective laser melting [22] and selective binding processes [23,24]. For structures built with these technologies, the thickness of the shell is often a simple result of the line width of the printer and is hence not specified as an infill volume. Here, the designer is interested in the resulting shell geometry actually manufactured by the printer which can be predicted with our method provided that the line width is known. Such knowledge is of particular interest for the modeling and simulation of 3D-printed structures as the simulation software requires accurate volume geometries. In simulation software, the substantial reduction of the memory requirements due to the use of an analytic function instead of a high-resolution surface mesh is a substantial advantage as simulation software is typically memory bound and preprocessing the geometry often takes more time than the actual simulation. This is actually the authors’ main motivation as we are concerned with the simulation of air flow through 3D-printed TPMS with our lattice Boltzmann solver VirtualFluids [25,26,27,28]. Surface meshes are particularly ill-suited for large-scale periodic structures such that the search for a more economical representation of such structures was inevitable.

We note here that the proposed method is markedly different from discretized numerical methods for solving the eikonal equation such as the fast marching method described by Sethian [8]. The fast marching method depends on numerical differentiation, which is done on a grid. Even though Sethian points out that the fast marching method can be improved by using finite differences of higher order, this does not necessarily enhance the convergence order of the method. In fact, the convergence order is not directly related to the quality of the normal derivative as is seen from the fact that our method always builds on exact analytic derivations. Due to this utilization of exact derivatives, no discretization and hence no grid is required in the proposed method. This does not mean that our method is generally more efficient than a fast marching method, as discretization on a grid has its advantages. The analytically derived offset function can become extremely complicated. The proposed method is hence not to be seen as a replacement for the fast marching method but rather, as an additional, very precise tool for cases where the original implicit surface is efficiently expressed as an analytic function.

When the thickness of the structure is not constant but varies in space, our method reduces to second-order accuracy, which is still much more accurate than required by any perceivable application in additive manufacturing at the time.

Finally, we note that the extended Richardson extrapolation is by no means restricted to the method derived in this paper. The extended Richardson extrapolation provides a systematic procedure to derive methods of arbitrary accuracy starting from a poor approximation, provided that this approximation fulfills the requirement that the error can be expanded in series.