# TPMS2STEP: error-controlled and $C^{2}$ continuity-preserving translation of
TPMS models to STEP files based on constrained-PIA

Yaonaiming Zhao  Qiang Zou<EMAIL_ADDRESS>Guoyue Luo  Jiayu Wu
Sifan Chen  Depeng Gao  Minghao Xuan  Fuyu Wang  State Key Laboratory of
CAD$\&$CG, Zhejiang University, Hangzhou, 310027, China  Aviation Industry
Shenyang Aircraft Design and Research Institute, Liaoning, 110035, China

###### Abstract

Triply periodic minimal surface (TPMS) is emerging as an important way of
designing microstructures. However, there has been limited use of commercial
CAD/CAM/CAE software packages for TPMS design and manufacturing. This is
mainly because TPMS is consistently described in the functional representation
(F-rep) format, while modern CAD/CAM/CAE tools are built upon the boundary
representation (B-rep) format. One possible solution to this gap is
translating TPMS to STEP, which is the standard data exchange format of
CAD/CAM/CAE. Following this direction, this paper proposes a new translation
method with error-controlling and $C^{2}$ continuity-preserving features. It
is based on an approximation error-driven TPMS sampling algorithm and a
constrained-PIA algorithm. The sampling algorithm controls the deviation
between the original and translated models. With it, an error bound of
$2\epsilon$ on the deviation can be ensured if two conditions called
$\epsilon$-density and $\epsilon$-approximation are satisfied. The
constrained-PIA algorithm enforces $C^{2}$ continuity constraints during TPMS
approximation, and meanwhile attaining high efficiency. A theoretical
convergence proof of this algorithm is also given. The effectiveness of the
translation method has been demonstrated by a series of examples and
comparisons.

###### keywords:

CAD/CAM/CAE integration; TPMS models; STEP files; Error control; Constrained
NURBS approximation; $C^{2}$ continuity

††journal: SPM

##  1 Introduction

The triply periodic minimal surface (TPMS) is a porous cellular-like surface
with zero mean curvature [1]. TPMS structures have advantages such as high
porosity, high surface-to-volume ratios, and super mechanical properties [2],
etc. With these features, as well as recent advancements in additive
manufacturing technologies [3], come increasing applications of TPMS
structures in fields like aerospace [4], biomedical [5], and energy [6].

The design and manufacturing of TPMS structures require tight integration of
computer-aided design, manufacturing, and engineering (CAD/CAM/CAE) tools.
However, TPMS is usually described in the functional representation (F-rep)
format [7], while modern CAD/CAM/CAE software packages are built upon the
boundary representation (B-rep) format [8]. There are two possible solutions
to this problem: (1) F-rep-based integration using specialized CAD/CAM/CAE
software [2, 9, 10, 7]; (2) STEP-based integration using existing commercial
CAD/CAM/CAE software [11, 12, 13]. (STEP stands for the Standard for The
Exchange of Product data, and it is a B-rep-based neutral file format
supporting the solid model exchanges among CAD/CAM/CAE tools [14].)

Developing brand-new, specialized CAD/CAM/CAE software that can run stably and
efficiently in real industrial practice requires huge investment and a long
time of fine-tuning and user training [15]. This work thus opts for the second
strategy, which allows TPMS structures to be directly designed and made in
existing CAD/CAM/CAE software. Nevertheless, this cannot happen without a tool
that is able to translate TPMS models to STEP files. Existing methods related
to this topic may be classified into two major categories. The first one is to
approximate TPMS with surface meshes [12, 11]. The second category converts
TPMS to Non-Uniform Rational B-Splines (NURBS) surfaces [13]. Compared to
meshes, NURBS surfaces offer better accuracy and compactness in representing
shapes [16]. However, existing TPMS-to-NURBS methods are consistently limited
to simple cases, without considering error control over the approximation
process and continuity preservation in the translated model [13]. Without
controlled errors and smoothness, the translated model cannot be used by
downstream applications with high fidelity.

This paper presents our attempts to attain error control and $C^{2}$
continuity in TPMS-to-STEP translation for the most commonly used Gyroid,
Diamond, and Schwarz_P TPMS structures. The proposed method can guarantee that
the deviation between the original structure and the translated solid model
will not exceed a given error threshold. That is, an error bound of
$2\epsilon$ on the deviation can be ensured if two conditions called
$\epsilon$-density and $\epsilon$-approximation are satisfied. Besides error
controlling, the proposed method can preserve the $C^{2}$ continuity of
translated models. $C^{2}$ continuity constraints are first derived and then
incorporated into the traditional progressive iterative approximation (PIA)
algorithm [17] to attain an extended version of PIA, called constrained-PIA or
simply CPIA. CPIA not only provides the $C^{2}$ continuity feature but also
has high computational efficiency. This is particularly useful for cases
requiring high translation accuracy, which implies a large number of sampling
points to be approximated. We also prove that CPIA has a guaranteed
convergence. Note that CPIA is general enough to handle a higher order of
continuity, but this work focuses primarily on $C^{2}$ continuity.

The following sections begin with a review of related work in Sec. 2, then the
proposed method’s details in Sec. 3, then the validation of the method using a
series of examples and comparisons in Sec. 4, and finally conclusions on the
method’s advantages and limitations in Sec. 5.

##  2 Related Work

This section briefly discusses the literature related to TPMS modeling (Sec.
2.1) and NURBS approximation (Sec. 2.2). The TPMS modeling methods are
classified into two categories: implicit modeling and explicit modeling. The
process of NURBS approximation involves four steps: parameterization, knot
placement, weights assignment, and control points computation. Existing
methods for each step are reviewed.

###  2.1 TPMS modeling

Implicit methods. Implicit models use mathematical functions to represent the
shape of TPMS structures. Various mathematical functions have been used for
this purpose, including trigonometric equation-based method [1, 18], signed
distance field-based method [9], and boundary curves-based method [10, 19],
among others. Despite wide applications, these methods have some serious
limitations. For example, because of the limited expressiveness of
trigonometric equations, the trigonometric equation-based method can only
express a small part of TPMS structures. The signed distance field-based
method is general, but when it is converted to triangular meshes, sharp
features are often lost and low-quality meshes are generated. What’s worse,
implicit TPMS models may still need to be explicitized (i.e., to meshes) to
make downstream applications like interaction, analysis, and fabrication
easier [12, 11]. Mesh-free interaction, analysis, and fabrication of TPMS
models are possible, and considerable progress has been made, e.g., [3, 20].
This direction of research is promising, but they are essentially developing
another collection of CAD/CAM/CAE tools that completely leave B-rep techniques
out, which is never easy and requires a huge investment and a long time of
fine-tuning and user training [15].

Explicit methods. In contrast to implicit methods, explicit methods can make
full use of existing CAD/CAM/CAE software packages since they are B-rep-based.
Over the past few decades, various B-rep schemes have been employed to model
TPMS, including mesh surface-based method [12, 21, 22, 23, 24], subdivision
surface-based method [25, 26, 27], and NURBS surface-based method [13]. Among
them, the NURBS-based methods offer more potential for interfacing with
existing CAD/CAM/CAE tools due to the fundamental role of NURBS
curves/surfaces in CAD/CAM/CAE [28, 29]. In particular, the method presented
by Flores et al. [13] decomposes the problem of approximating TPMS models with
NURBS surfaces into several sub-problems of approximating a portion of TPMS
models and then assembling them to attain the overall approximation. This way
of working is easier than the original problem. Unfortunately, only simple
cases were considered, and how to control approximation error and preserve
continuity in assembling individual approximations is not quite known. The
proposed method also employs this decomposition-then-recombination framework
but with new error-controlling and $C^{2}$ continuity-preserving features.

###  2.2 NURBS approximation

The underlying technique of TPMS-to-STEP translation is NURBS surface
approximation, which is in turn related to parameterization, knot placement,
weight assignment, and computation of control points. For this reason, this
work also gives a brief review of each of these procedures, as follows. A
thorough review of them can be found in existing literature, e.g., [30].

Parameterization. Parametrization means associating each data point with a
parameter value. Heuristic methods are the dominant means of parameterization,
e.g., uniform [30], chord length [31], centripetal [32], universal [33], and
Foley-Nielson [34], to name a few. Such heuristic methods have also been used
as an initial guess for further parameterization, e.g., [35].

Knot placement. Knot placement refers to the choice of the number of knots and
their locations for constructing B-spline bases. Like parameterization,
heuristic methods like the averaging technique (AVG) and the knot placement
technique (KTP) are often the method of choice in practice [30, 36]. These
methods have also been used as an initial guess and then optimized through
iterative knot insertion [37, 38].

Weight assignment. Weights are used to attain more expressibility in
representing shapes with NURBS. Weights are usually added in an ad-hoc manner,
e.g., the uniform method [30] and the curvature-based method [39].

The parameterization, knot placement, and weight assignment methods mentioned
above provide a variety of options in NURBS approximation. In this paper, to
simplify the NURBS approximation problem in TPMS2STEP and make it easier to
implement, uniform parameterization, uniform knot placement, and uniform
weight assignment methods are chosen.

Control points computation. The control points computation is typically
expressed in terms of a least square optimization problem, which boils down to
solving a linear system [13]. This method is time-consuming and becomes
infeasible when the number of the data points to be approximated is very large
[17]. To address this issue, iterative methods have been used instead, such as
those presented in [40, 41]. Among them, the progressive iterative
approximation (PIA) method [17] may be the easiest to implement and provides
more flexibility to incorporate geometric constraints. Thus, it will be used
as the backbone algorithm for developing the error-controlled and $C^{2}$
continuity-preserving TPMS approximation algorithm in this paper.

As can be seen from the above review, NURBS approximation is a well-
established field, as long as the data points are ready and of good quality.
In this regard, this work does not intend to add new results to NURBS
approximation but to adapt the existing methods stated above to the special
needs of TPMS-to-STEP translation. That is, we focus primarily on two
problems: (1) how to prepare quality data points for NURBS approximation so
that the final results have controlled errors; (2) how to ensure that the
continuity of original TPMS models is not broken in the approximation.

##  3 Methods

###  3.1 TPMS and NURBS preliminaries

To better present TPMS2STEP, some basic knowledge of TPMS and NURBS, as well
as their notations, is first described. TPMS refers to a surface with zero
mean curvature arranged periodically in 3D space. To design a TPMS solid model
with non-zero thickness (Fig. 1), the equations for offset primitive surfaces
can be derived from the Enneper–Weierstrass representation [42, 43, 44] :

$\displaystyle\begin{aligned}
&x=e^{i\theta}Re\left[\int_{\omega_{0}}^{\omega}(1-\tau^{2})R(\tau)d\tau\right]+\frac{n_{x}d}{\sqrt[]{n_{x}^{2}+n_{y}^{2}+n_{z}^{2}}}\\\
&y=e^{i\theta}Re\left[\int_{\omega_{0}}^{\omega}i(1+\tau^{2})R(\tau)d\tau\right]+\frac{n_{y}d}{\sqrt[]{n_{x}^{2}+n_{y}^{2}+n_{z}^{2}}}\\\
&z=e^{i\theta}Re\left[\int_{\omega_{0}}^{\omega}2\tau
R(\tau)d\tau\right]+\frac{n_{z}d}{\sqrt[]{n_{x}^{2}+n_{y}^{2}+n_{z}^{2}}}\end{aligned}$
(1)

where $i$ is an imaginary unit. $Re$ extracts the real part of a complex
number. $d$ is the offset value. $\theta$ is the Bonnet angle. $R(\tau)$ is
the Weierstrass function with a complex variant $\tau$. $\omega_{0}$ and
$\omega$ are the two end points of the integral path in Enneper–Weierstrass
representation. The normal vectors $n=(n_{x},n_{y},n_{z})$ are defined as
[45]:

$\displaystyle\begin{aligned} &n_{x}=2Im(\phi_{2}\overline{\phi_{3}})\\\
&n_{y}=2Im(\phi_{3}\overline{\phi_{1}})\\\
&n_{z}=2Im(\phi_{1}\overline{\phi_{2}})\end{aligned}$

where $\phi_{1}$, $\phi_{2}$, and $\phi_{3}$ are the first-order derivatives
of Enneper–Weierstrass representation. $\overline{\phi_{1}}$,
$\overline{\phi_{2}}$, and $\overline{\phi_{3}}$ are the conjugation of
$\phi_{1}$, $\phi_{2}$, and $\phi_{3}$ respectively. $Im$ refers to the
operator taking the imaginary part of a complex number.

Figure 1: Illustration of TPMS solid models for Gyroid (left), Diamond
(middle), and Schwarz_P (right).

A NURBS surface is described by:

$\mathbf{S}(u,v)=\frac{\sum_{i=0}^{n}\sum_{j=0}^{m}\mathbf{B}_{i,p}(u)\mathbf{B}_{j,q}(v)w_{i,j}\mathbf{P}_{i,j}}{\sum_{i=0}^{n}\sum_{j=0}^{m}\mathbf{B}_{i,p}(u)\mathbf{B}_{j,q}(v)w_{i,j}}$
(2)

where $u$ and $v$ are parameters in $[0,1]$.
$\mathbf{B}_{i,p}(u),i=0,\cdots,n,$ and $\mathbf{B}_{j,q}(v),j=0,\cdots,m,$
are the B-spline basis functions of degree p and q.
$\mathbf{P}_{i,j},i=0,\cdots,n,j=0,\cdots,m$ are control points. $w_{i,j}$
means the weight of the corresponding control point $\mathbf{P}_{i,j}$. The
basis functions $\mathbf{B}_{i,p}(u)$ and $\mathbf{B}_{j,q}(v)$ are generated
using a knot vector $U=\\{0=u_{0}<u_{1}<\cdots<u_{n}=1\\}$ in $u$ direction
and a knot vector $V=\\{0=v_{0}<v_{1}<\cdots<v_{m}=1\\}$ in $v$ direction, as
follows:

$\displaystyle\mathbf{B}_{i,0}(u)=\left\\{\begin{array}[]{rcl}1&if\;u_{i}\leq
u\leq u_{i+1}\\\ 0&otherwise\end{array}\right.$
$\displaystyle\mathbf{B}_{i,p}(u)=\frac{u-u_{i}}{u_{i+p}-u_{i}}\mathbf{B}_{i,p-1}(u)+\frac{u_{i+p+1}-u}{u_{i+p+1}-u_{i+1}}\mathbf{B}_{i+1,p-1}(u)$

###  3.2 The TPMS2STEP pipeline

The proposed TPMS-to-STEP translation pipeline is shown in Fig. 2. It begins
with four input parameters: the type of TPMS, the offset value for turning a
TPMS surface into a solid, the approximation tolerance, and possibly the
scaling factor varying with location. This input can be viewed as a parametric
representation of the intended TPMS structure. This parametric representation
is converted to NURBS surfaces with regard to the approximation tolerance
specified by the user and the $C^{2}$ continuous constraints intrinsic to TPMS
itself. This is where the proposed CPIA comes into play. Because TPMS is
highly periodic, there is no need for approximating all surfaces bounding the
TPMS structure but only several representative surfaces (to be called
primitive surfaces in what follows). Having primitive surfaces in place, we
create several copies of them and reposition them to the proper locations,
followed by a sewing step to form a watertight, $C^{2}$ continuous B-rep
model. The last step is to apply the varying scaling factor, if any, to this
B-rep model and then write the scaled model to a .step file according to the
STEP format.

Figure 2: The proposed TPMS-to-STEP translation pipeline.

In the above process of constructing NURBS surfaces, the error between the
original surface and the approximation surface needs to be controlled so that
the resulting STEP file can be used by downstream applications with high
fidelity. Because it is hard, if not impossible, to simultaneously carry out
TPMS sampling and NURBS fitting with respect to an error limit, a virtual
intermediate surface is created to divide the overall error into two decoupled
errors: one from the TPMS samples to the virtual intermediate surface, and the
other from the virtual intermediate surface to the approximate NURBS surface.
Such a decoupling of errors allows partial error control. That is, if the
error from the original TPMS model to this intermediate surface is $\epsilon$,
and the other error from the final approximation surface to this intermediate
surface is $\epsilon$, then the overall error must not exceed $2\epsilon$, as
shown in Fig. 3. Then this error bound $2\epsilon$ can be forced to equal the
user-specified tolerance, thereby simplifying the error-controlling task. For
each partial error control, an error estimation model is to be developed so
that the sampling can be not only driven by the user-specified tolerance but
also done adaptively.

Figure 3: 2D illustration of the virtual intermediate surface method. The
original surface and approximation surfaces are constrained within the error
band (the dashed lines), and therefore a $2\epsilon$ error bound on their
deviation.

The other major problem in NURBS approximation is $C^{2}$ continuity. As
already noted, the PIA algorithm [17] is to be used to carry out the
approximation. However, a direct application of PIA does not help here because
there are complex geometric constraints involved and they can change the
original PIA iteration scheme into a new form. A modified version of PIA,
i.e., CPIA, is proposed for this reason. In Sec. 3.4, we first derive those
geometric constraints regarding $C^{2}$ continuity and then incorporate them
into the PIA algorithm to attain $C^{2}$ continuous TPMS approximation. A
theoretical convergence proof for CPIA’s new iteration scheme is given in Sec.
3.5.

###  3.3 Error-controlled TPMS sampling

As already noted, using a virtual intermediate surface, an error bound of
$2\epsilon$ can be attained on the deviation between the original surface and
the approximation surface. The essential tasks are then constructing the
virtual intermediate surface and controlling the errors from the original
surface to the virtual intermediate surface (abbr. as O-to-V) and from the
approximation surface to the virtual intermediate surface (abbr. as A-to-V) so
that they both are under $\epsilon$. In this work, triangular meshes are
chosen to construct virtual intermediate surfaces. Actually, constructing the
virtual intermediate surface and controlling O-to-V error can be done in a
combined way. This is because the O-to-V error is equivalent to the V-to-O
error. We only need to sample the original surface in such a density that the
mesh from these sample points has a distance of $\epsilon$ to the original
surface. We call this $\epsilon$-density. Similarly, the approximation that
gives an A-to-V error of $\epsilon$ is to be called $\epsilon$-approximation.

Í

####  3.3.1 $\epsilon$-density

As can be seen from the zoomed-in illustration of Fig. 4, $\epsilon$-density
is to determine the $\Delta$ such that the local shape of the orginal surface
$S_{O}(u,v)$ is bounded by the polyhedron generated through offsetting the
triangles $(v_{1},v_{2},v_{3})$ and $(v_{1},v_{3},v_{4})$ by $\epsilon$, where
$v_{1}=S_{O}(u_{0},v_{0})$, $v_{2}=S_{O}(u_{0}+\Delta,v_{0})$,
$v_{3}=S_{O}(u_{0},v_{0}+\Delta)$, and
$v_{4}=S_{O}(u_{0}+\Delta,v_{0}+\Delta)$.

The density parameter $\Delta$ and the offset value $\epsilon$ are related by
the following equation (the Fillip’s algorithm [46]):

$\epsilon=\frac{(M_{1}+2M_{2}+M_{3})}{8}\Delta^{2}$ (3)

where $M_{1}$, $M_{2}$, and $M_{3}$ are the maximum second-order partial
derivatives of $S_{O}(u,v)$, as follows:

$\displaystyle M_{1}=\max_{(u,v)\in D}(\left|\frac{\partial^{2}x}{\partial
u^{2}}\right|,\left|\frac{\partial^{2}y}{\partial
u^{2}}\right|,\left|\frac{\partial^{2}z}{\partial u^{2}}\right|)$
$\displaystyle M_{2}=\max_{(u,v)\in D}(\left|\frac{\partial^{2}x}{\partial
u\partial v}\right|,\left|\frac{\partial^{2}y}{\partial u\partial
v}\right|,\left|\frac{\partial^{2}z}{\partial u\partial v}\right|)$
$\displaystyle M_{3}=\max_{(u,v)\in D}(\left|\frac{\partial^{2}x}{\partial
v^{2}}\right|,\left|\frac{\partial^{2}y}{\partial
v^{2}}\right|,\left|\frac{\partial^{2}z}{\partial v^{2}}\right|)$

$D$ is the parameter domain of interest, and the method to estimate the
maximum second-order partial derivatives is given in Section 3 of the
supplementary material.

Based on Eq. (3), the $\epsilon$-density algorithm consists of three major
steps:

  1. 1.

Substitute the given error threshold $\epsilon$ to get the sampling density
parameter $\Delta$;

  2. 2.

Discretize the parameter domain of $S_{O}(u,v)$ into a grid with size
$\Delta$; and

  3. 3.

Map the parametric grid to 3D space through $S_{O}(u,v)$ and then mesh it into
the intended virtual intermediate surface.

Figure 4: 2D illustration of $\epsilon$-density.

####  3.3.2 $\epsilon$-approximation

Because it is hard to estimate the error from the approximation surface to the
virtual intermediate surface before approximation, the error is controlled
with an indirect approach. That is, the error is estimated after approximation
and is utilized to refine the sampling. This refinement increases the sampling
density if the error is beyond $\epsilon$. For error estimation, the algorithm
proposed by Zheng et al. [47] is adopted due to its tightness. The error
estimation equation is:

$\displaystyle\epsilon_{ij}$
$\displaystyle=\frac{L_{grid}^{2}n^{2}(\mathbf{A}_{ij}+r\mathbf{G}_{ij})}{4\mathrm{inf}\\{w(u,v)\\}+L_{grid}^{2}n^{2}\mathbf{G}_{ij}},0\leq
i\leq n-1,0\leq j\leq m-1$ (4) $\displaystyle\mathbf{A}_{ij}$
$\displaystyle=\parallel
w_{i,j}\mathbf{P}_{i,j}-w_{i+1,j}\mathbf{P}_{i+1,j}-w_{i,j+1}\mathbf{P}_{i,j+1}+w_{i+1,j+1}\mathbf{P}_{i+1,j+1}\parallel$
$\displaystyle\mathbf{G_{ij}}$ $\displaystyle=\mid
w_{i+1,j+1}-w_{i+1,j}-w_{i,j+1}+w_{i,j}\mid$

where $\epsilon_{ij}$ is the estimated tolerance in the $(i,j)$th grid of
control points, $L_{grid}$ the distance between each two grid points, $m$ and
$n$ the number of control points in $u$ and $v$ directions plus one, and
$r=\max_{0\leq i\leq n\atop 0\leq j\leq m}\parallel\mathbf{P}_{ij}\parallel$.

After error estimation, the adaptive refinement of sampling density is carried
out by adding a row and a column of sample points locally in the grids where
the local error is larger than $\epsilon$. The above refinement is repeated
until the approximation error is below $\epsilon$.

###  3.4 $C^{2}$ continuity-preserving TPMS approximation

Approximating data points ${\mathbf{Q}_{ij}}(i=0,\dots,m_{1},j=0,\dots,m_{2})$
with a NURBS surface is to solve the following optimization problem [13]:

$\min_{{\mathbf{P}_{00},\cdots,\mathbf{P}_{n_{1}n_{2}}}}\sum_{i=0}^{m1}\sum_{j=0}^{m2}||\mathbf{Q}_{ij}-\sum_{k=0}^{n1}\sum_{l=0}^{n2}\mathbf{B}(u_{i})\mathbf{B}(v_{j})w_{kl}\mathbf{P}_{kl}||^{2}_{2}$
(5)

where $u_{i}$ and $v_{j}$ are the parameter values corresponding to the data
point $\mathbf{Q}_{ij}$. The parameter values are generated using the
following procedures. First, the fundamental patch is sampled into point
clouds, both in the 2D parameter domain and the 3D space. Then the parameter
domain of the fundamental patch is mapped to a triangular region. This region
is extended to a rectangular region by the symmetry property, and the same for
2D parametric samples, which results in a quasi-grid of parametric samples.
Finally, those parametric samples are associated with the regular grid points
in [0,1]×[0,1], correspondingly.

The above optimization problem is a least square problem, and its solution can
be obtained by solving the following linear system:

$(\mathbf{B}^{T}\mathbf{B})\mathbf{P}=\mathbf{B}^{T}\mathbf{Q}$ (6)

where $\mathbf{P}=[\mathbf{P}_{00},\cdots,\mathbf{P}_{n_{1}n_{2}}]^{T}$,
$\mathbf{Q}=[\mathbf{Q}_{00},\cdots,\mathbf{Q}_{m_{1}m_{2}}]^{T}$, and

$\displaystyle\mathbf{B}=\begin{bmatrix}\mathbf{B}_{0,p}(u_{0})\mathbf{B}_{0,p}(v_{0})&\dots&\mathbf{B}_{n_{1},p}(u_{0})\mathbf{B}_{n_{2},p}(v_{0})\\\
\vdots&\ddots&\vdots\\\
\mathbf{B}_{0,p}(u_{m_{1}})\mathbf{B}_{0,p}(v_{m_{2}})&\dots&\mathbf{B}_{n_{1},p}(u_{m_{1}})\mathbf{B}_{n_{2},p}(v_{m_{2}})\end{bmatrix}$

If a high TPMS translation accuracy is required (i.e., large-scale sampling
points to be approximated), the matrix $\mathbf{B}$ in Eq. (6) is large, and
solving this equation is time-consuming. For this reason, the fast PIA
algorithm is used instead of solving the equation directly. Its basic steps
are summarized below.

The initial step of PIA is constructing a NURBS surface using the data points
$\mathbf{Q}$ directly as control points, as follows:

$\mathbf{P}^{0}(u,v)=\sum_{i=0}^{m_{1}}\sum_{j=0}^{m_{2}}\mathbf{B}_{i}(u)\mathbf{B}_{j}(v)w_{ij}\mathbf{P}^{0}_{ij},0\leq
u,v\leq 1$ (7)

where $\mathbf{P}^{0}_{ij}=\mathbf{Q}_{ij}$. The second step is computing the
deviation from each ${\mathbf{Q}}_{ij}$ to its associative point on surface
$\mathbf{P}^{0}(u,v)$, as follows:

$\mathbf{\Delta}^{0}_{ij}=\mathbf{Q}_{ij}-\mathbf{P}^{0}(u_{i},v_{j})$ (8)

The third step is using this deviation to update the initial control points in
Eq. (7):

$\mathbf{P}^{1}_{ij}=\mathbf{P}^{0}_{ij}+\mathbf{\Delta}^{0}_{ij}$ (9)

Repeating Eqs. (8) and (9) yields a sequence of surfaces
$\\{\mathbf{P}^{0}(u,v),\mathbf{P}^{1}(u,v),\cdots,\mathbf{P}^{k}(u,v)\\}$.
The limit of this surface sequence is equivalent to solving Eq. (6), as proved
by Lin et al. [17]. It is easy to see that this iterative scheme is highly
parallel, thereby an efficient method.

A direct application of PIA is not enough for TPMS approximation. $C^{2}$
continuity constraints for smoothly sewing primitive surfaces need to be
derived and then incorporated into the PIA algorithm, leading to a new
iteration scheme (i.e., CPIA). In the following, we first derive those
constraints and then present CPIA’s details, followed by the convergence proof
of its iteration scheme.

####  3.4.1 $C^{2}$ continuity constraints

To preserve $C^{2}$ continuity in approximation, geometric constraints need to
be added where surfaces meet. To find the junctions of surfaces, we carefully
considered the surface assembly methods used in this paper (i.e., the methods
in [13, 42, 44]). As is shown in Fig. 5, the assembled offset surfaces are
constructed by the copy and rigid transformation of two primitive surfaces
with $offset=+\delta$ and $offset=-\delta$. The specific rigid transformation
matrices for assembly can be found in Eq. (6), Eq. (16), and Eq. (19) of [13]
for Gyroid, Table 1 in [42] for Diamond, and Table 1 in [44] for Schwarz_P.
(We omit the details of the matrices here because they are not new results.
Also, the specific entries in those matrices are not important, nor do they
affect the presentation of this work.) Because these surfaces are arranged
regularly, the constraints at all the junctions of surfaces can be reduced to
four fundamental situations shown in Table 1, with the number of edges denoted
as $e_{1}^{+},e_{2}^{+},e_{3}^{+},e_{4}^{+}$ and
$e_{1}^{-},e_{2}^{-},e_{3}^{-},e_{4}^{-}$. Other situations can be obtained by
rigid transformation of these four situations.

Figure 5: Illustration of the $C^{2}$ continuity-preserving method: (a) Gyroid; (b) Diamond; (c) Schwarz_P. The edges of the surfaces are divided into two, two, and four types with different colors for Gyroid, Diamond, and Schwarz_P TPMS respectively. Only the edges with the same color can be spliced together during assembly. Table 1: Constraints of the edges of primitive surfaces for Gyroid, Diamond, and Schwarz_P. $\mathbf{T}_{ig}(i=1,2)$, $\mathbf{T}_{id}(i=1,2)$, and $\mathbf{T}_{ip}(i=1,2,3,4)$ are given in Section 1 of the supplementary material. Model Type | Constraints
---|---
Gyroid | $e_{1}^{+}=\mathbf{T}_{1g}e_{2}^{-}$ $e_{2}^{+}=\mathbf{T}_{1g}^{-1}e_{1}^{-}$ $e_{3}^{+}=\mathbf{T}_{2g}e_{4}^{-}$ $e_{4}^{+}=\mathbf{T}_{2g}^{-1}e_{3}^{-}$
Diamond | $e_{1}^{+}=\mathbf{T}_{1d}e_{4}^{-}$ $e_{2}^{+}=\mathbf{T}_{2d}e_{3}^{-}$ $e_{3}^{+}=\mathbf{T}_{2d}^{-1}e_{2}^{-}$ $e_{4}^{+}=\mathbf{T}_{1d}^{-1}e_{1}^{-}$
Schwarz_P | $e_{1}^{+}=\mathbf{T}_{1p}e_{1}^{+}$ $e_{2}^{+}=\mathbf{T}_{2p}e_{2}^{+}$ $e_{3}^{+}=\mathbf{T}_{3p}e_{3}^{+}$ $e_{4}^{+}=\mathbf{T}_{4p}e_{4}^{+}$

For each constraint in Table 1, a constraint equation is derived to adjust the
corresponding control points of the approximation surface to keep the first-
order and second-order derivatives at the junctions to be the same [48].
Because the edge-edge connection relationships for the three types of TPMS are
different, their constraints are different and thus considered separately.

Gyroid. For Gyroid TPMS, the two primitive surfaces used to assemble the
offset surface have the same shape (shown in Fig. 5a). That is, the primitive
surface with $offset=-\delta$ can be obtained by taking a specific rigid
transformation, denoted as $\mathbf{T}_{3g}$, to the primitive surface with
$offset=+\delta$. $\mathbf{T}_{3g}$ is given in Section 1 of the supplementary
material. Based on this observation, the splicing of two primitive surfaces
can be considered as splicing one primitive surface with another primitive
surface copied and transformed by itself, and the constraints for different
edges of the single primitive surface are derived from the constraints for
edges of two primitive surfaces in Table 1, given as:

$e_{1}^{+}=\mathbf{T}_{4g}e_{1}^{+},\quad
e_{2}^{+}=\mathbf{T}_{5g}e_{2}^{+},\quad
e_{3}^{+}=\mathbf{T}_{6g}e_{3}^{+},\quad e_{4}^{+}=\mathbf{T}_{7g}e_{4}^{+}$
(10)

where $\mathbf{T}_{4g}$-$\mathbf{T}_{7g}$ are derived from
$\mathbf{T}_{1g}$-$\mathbf{T}_{3g}$ and they are given in Section 1 of the
supplementary material.

To ensure $C^{2}$ continuity at the junctions, part of the control points of
the approximation surface that control the first-order and second-order
derivatives of the two connected edges are extracted and adjusted for
consistency. Specifically, for the four constraints in Eq. 10, if the control
points of the approximation surface are denoted as $\mathbf{P}$ and the group
of matrices used to extract the corresponding control points are denoted as
$\mathbf{M}_{ijg}(i=1,2,3,4,5$ and $j=1,2,3,4)$, then the constraint equations
are derived as:

$\displaystyle\begin{aligned}
&(\mathbf{M}_{2ig}-\mathbf{M}_{1ig})\mathbf{P}\mathbf{T}_{(i+3)g}=(\mathbf{M}_{2ig}-\mathbf{M}_{1ig})\mathbf{P}\\\
&(3\mathbf{M}_{4ig}-2\mathbf{M}_{3ig}-\mathbf{M}_{5ig})\mathbf{P}\mathbf{T}_{(i+3)g}=(3\mathbf{M}_{4ig}-2\mathbf{M}_{3ig}-\mathbf{M}_{5ig})\
\mathbf{P}\\\ &i=1,2,3,4\end{aligned}$ (11)

where the left side of the equation indicates the first-order and second-order
derivatives of one edge and the right side of the equation indicates the
first-order or second-order derivatives of the other edge, and the same
coefficients on both sides of the equations to calculate the first-order or
second-order derivatives of NURBS surfaces are eliminated to simplify the
form. $\mathbf{M}_{ijg}$ $(i=1,2,3,4,5$ and $j=1,2,3,4)$ are used to extract
the first, second, and third rows of control points at the four edges of the
approximation surface. Specifically,
$(\mathbf{M}_{2ig}-\mathbf{M}_{1ig})\mathbf{P}$ $(i=1,2,3,4)$ are utilized to
calculate the first-order derivative of the surface, and
$(3\mathbf{M}_{4ig}-2\mathbf{M}_{3ig}-\mathbf{M}_{5ig})\mathbf{P}$
$(i=1,2,3,4)$ are utilized to calculate the second-order derivative of the
surface. Their details are provided in Section 1 of the supplementary
material.

Diamond.

For Diamond TPMS, the two primitive surfaces for assembling have different
shapes (shown in Fig. 5b). To ensure $C^{2}$ continuity with the same approach
for Gyroid, the continuity for each constraint (shown in Table 1) is
considered during the approximation of these two primitive surfaces. Similar
to Gyroid, the $C^{2}$ continuity constraint equations are derived as:

$\displaystyle\begin{aligned}
&(\mathbf{M}_{24d}-\mathbf{M}_{14d})\mathbf{N}_{i}\mathbf{P}\mathbf{T}_{1d}=(\mathbf{M}_{21d}-\mathbf{M}_{11d})\mathbf{N}_{3-i}\mathbf{P}\\\
&(3\mathbf{M}_{44d}-2\mathbf{M}_{34d}-\mathbf{M}_{54d})\mathbf{N}_{i}\mathbf{P}\mathbf{T}_{1d}=(3\mathbf{M}_{41d}-2\mathbf{M}_{31d}-\mathbf{M}_{51d})\mathbf{N}_{3-i}\mathbf{P}\\\
&(\mathbf{M}_{23d}-\mathbf{M}_{13d})\mathbf{N}_{i}\mathbf{P}\mathbf{T}_{2d}=(\mathbf{M}_{22d}-\mathbf{M}_{12d})\mathbf{N}_{3-i}\mathbf{P}\\\
&(3\mathbf{M}_{43d}-2\mathbf{M}_{33d}-\mathbf{M}_{53d})\mathbf{N}_{i}\mathbf{P}\mathbf{T}_{2d}=(3\mathbf{M}_{42d}-2\mathbf{M}_{32d}-\mathbf{M}_{52d})\mathbf{N}_{3-i}\mathbf{P}\\\
&(i=1,2)\end{aligned}$ (12)

where $\mathbf{P}=\begin{pmatrix}\mathbf{P}_{1}\\\
\mathbf{P}_{2}\end{pmatrix}$. $\mathbf{P}_{1}$ and $\mathbf{P}_{2}$ are the
homogeneous coordinates of control points for the two approximation surfaces.
$\mathbf{M}_{ijd}$ $(i=1,2,3,4,5$ and $j=1,2,3,4)$ have the same function as
$\mathbf{M}_{ijg}$ $(i=1,2,3,4,5$ and $j=1,2,3,4)$ for Gyroid.
$\mathbf{N}_{1}$ and $\mathbf{N}_{2}$ are used to assemble the control points
$\mathbf{P}_{1}$ and $\mathbf{P}_{2}$ to $\mathbf{P}$. Specifically,
$\mathbf{N}_{1}\mathbf{P}=\begin{pmatrix}\mathbf{P}_{1}\\\
\mathbf{0}\end{pmatrix}$ and
$\mathbf{N}_{2}\mathbf{P}=\begin{pmatrix}\mathbf{P}_{2}\\\
\mathbf{0}\end{pmatrix}$. Matrices $\mathbf{T}_{id}$ $(i=1,2)$ are the rigid
transformation matrices given in Table 1. Their details can be found in
Section 1 of the supplementary material.

Schwarz_P. For Schwarz_P TPMS, only one primitive surface is needed to
assemble the offset surface (shown in Fig. 5c). That is, the offset surface
with $offset=-\delta$ can be obtained by copying and transforming the
primitive surface with $offset=+\delta$. The same approach to ensure $C^{2}$
continuity for Gyroid is used here for each of the constraints in Table 1 and
the constraint equations are derived as:

$\displaystyle\begin{aligned}
&(\mathbf{M}_{2ip}-\mathbf{M}_{1ip})\mathbf{P}\mathbf{T}_{ip}=(\mathbf{M}_{2ip}-\mathbf{M}_{1ip})\mathbf{P}\\\
&(3\mathbf{M}_{4ip}-2\mathbf{M}_{3ip}-\mathbf{M}_{5ip})\mathbf{P}\mathbf{T}_{ip}=(3\mathbf{M}_{4ip}-2\mathbf{M}_{3ip}-\mathbf{M}_{5ip})\mathbf{P}\\\
&i=1,2,3,4\end{aligned}$ (13)

where $\mathbf{P}$ is the vector of control points. $\mathbf{M}_{ijp}$
$(i=1,2,3,4,5$ and $j=1,2,3,4)$ are the same as $\mathbf{M}_{ijd}$
$(i=1,2,3,4,5$ and $j=1,2,3,4)$ for Diamond. Their details have been put in
Section 1 of the supplementary material.

For $C^{2}$ continuity-preserving, a general problem is that when there are
extraordinary points, it is hard to achieve high-order parametric continuity.
However, the three types of TPMS structures considered in this paper are
intrinsic $C^{2}$ continuous, effectively avoiding the issue of extraordinary
points.

####  3.4.2 The CPIA method

CPIA has the same major steps as PIA, but with a different iterative scheme.
The first step of CPIA is constructing a NURBS surface utilizing the sample
points $\mathbf{Q}$ as control points, denoted as
$\mathbf{P}^{0}(u_{i},v_{j})$. The second step is computing the update term
$\Delta$ for $\mathbf{P}^{0}(u_{i},v_{j})$. In PIA, this term refers to the
deviation from each sample point ${\mathbf{Q}}_{ij}$ to its associative point
on the approximation surface $\mathbf{P}^{0}(u,v)$. However, in CPIA, an
additional constraint term should also be considered to ensure that the
approximation surface satisfies the continuity constraints in Eq. (11), Eq.
(12), and Eq. (13).

For Gyroid, this term is derived from Eq. (11) by calculating the difference
between the two sides of the equation. This difference is actually the
deviation between the theoretical value that meets the constraints and the
actual value. Then a factor $\frac{1}{2}$ is added to distribute the deviation
into two equal parts since the deviation is caused by two sides of the
surfaces at their junction. The difference is added to the original PIA update
term (i.e., Eq. (8)) to obtain the new CPIA update term:

$\displaystyle\begin{aligned}
\mathbf{\Delta}^{0}_{ij}&=\mathbf{Q}_{ij}-\mathbf{P}^{0}(u_{i},v_{j})\\\
&+\frac{1}{2}\sum_{i=1}^{4}((\mathbf{M}_{2ig}-\mathbf{M}_{1ig})\mathbf{P}^{0}(u_{i},v_{j})\mathbf{T}_{(i+3)g}-(\mathbf{M}_{2ig}-\mathbf{M}_{1ig})\mathbf{P}^{0}(u_{i},v_{j}))\\\
&+\frac{1}{2}\sum_{i=1}^{4}((3\mathbf{M}_{4ig}-2\mathbf{M}_{3ig}-\mathbf{M}_{5ig})\mathbf{P}^{0}(u_{i},v_{j})\mathbf{T}_{(i+3)g}\\\
&-(3\mathbf{M}_{4ig}-2\mathbf{M}_{3ig}-\mathbf{M}_{5ig})\
\mathbf{P}^{0}(u_{i},v_{j}))\end{aligned}$ (14)

For Diamond, the additional constraint term is derived from Eq. (12) with
similar reasoning: calculating the difference between the two sides of the
equation and multiplying it with a factor $\frac{1}{2}$, and then adding it to
Eq. (8), yielding:

$\displaystyle\begin{aligned}
\mathbf{\Delta}^{0}_{ij}&=\mathbf{Q}_{ij}-\mathbf{P}^{0}(u_{i},v_{j})\\\
&+\frac{1}{2}\sum_{i=1}^{2}((\mathbf{M}_{24d}-\mathbf{M}_{14d})\mathbf{N}_{i}\mathbf{P}^{0}(u_{i},v_{j})\mathbf{T}_{1d}-(\mathbf{M}_{21d}-\mathbf{M}_{11d})\mathbf{N}_{3-i}\mathbf{P}^{0}(u_{i},v_{j}))\\\
&+\frac{1}{2}\sum_{i=1}^{2}((3\mathbf{M}_{44d}-2\mathbf{M}_{34d}-\mathbf{M}_{54d})\mathbf{N}_{i}\mathbf{P}^{0}(u_{i},v_{j})\mathbf{T}_{1d}\\\
&-(3\mathbf{M}_{41d}-2\mathbf{M}_{31d}-\mathbf{M}_{51d})\mathbf{N}_{3-i}\mathbf{P}^{0}(u_{i},v_{j}))\\\
&+\frac{1}{2}\sum_{i=1}^{2}((\mathbf{M}_{23d}-\mathbf{M}_{13d})\mathbf{N}_{i}\mathbf{P}^{0}(u_{i},v_{j})\mathbf{T}_{2d}-(\mathbf{M}_{22d}-\mathbf{M}_{12d})\mathbf{N}_{3-i}\mathbf{P}^{0}(u_{i},v_{j}))\\\
&+\frac{1}{2}\sum_{i=1}^{2}((3\mathbf{M}_{43d}-2\mathbf{M}_{33d}-\mathbf{M}_{53d})\mathbf{N}_{i}\mathbf{P}^{0}(u_{i},v_{j})\mathbf{T}_{2d}\\\
&-(3\mathbf{M}_{42d}-2\mathbf{M}_{32d}-\mathbf{M}_{52d})\mathbf{N}_{3-i}\mathbf{P}^{0}(u_{i},v_{j}))\end{aligned}$
(15)

For Schwarz_P, the update term follows directly as:

$\displaystyle\begin{aligned}
\mathbf{\Delta}^{0}_{ij}&=\mathbf{Q}_{ij}-\mathbf{P}^{0}(u_{i},v_{j})\\\
&+\frac{1}{2}\sum_{i=1}^{4}((\mathbf{M}_{2ip}-\mathbf{M}_{1ip})\mathbf{P}^{0}(u_{i},v_{j})\mathbf{T}_{ip}-(\mathbf{M}_{2ip}-\mathbf{M}_{1ip})\mathbf{P}^{0}(u_{i},v_{j}))\\\
&+\frac{1}{2}\sum_{i=1}^{4}((3\mathbf{M}_{4ip}-2\mathbf{M}_{3ip}-\mathbf{M}_{5ip})\mathbf{P}^{0}(u_{i},v_{j})\mathbf{T}_{ip}\\\
&-(3\mathbf{M}_{4ip}-2\mathbf{M}_{3ip}-\mathbf{M}_{5ip})\mathbf{P}^{0}(u_{i},v_{j}))\end{aligned}$
(16)

Having obtained the update terms, the third step utilizes them to update the
control points:

$\mathbf{P}^{1}_{ij}=\mathbf{P}^{0}_{ij}+\mathbf{\Delta}^{0}_{ij}$ (17)

Eqs. (14)-(17) will be repeated until the approximation error is below a given
tolerance. Algorithm 1 gives more details on these repetitive procedures.

Algorithm 1 The CPIA Algorithm

1:An ordered point set $\\{\mathbf{Q}_{ij}\\}^{m_{1},m_{2}}_{i=0,j=0}$ and
corresponding parameters $\\{\xi_{i},\zeta_{j}\\}^{m_{1},m_{2}}_{i=0,j=0}$,
weight of control points $w_{ij}$, knot vectors
$U=\\{0=u_{0}<u_{1}<\cdots<u_{n}=1\\}$ and
$V=\\{0=v_{0}<v_{1}<\cdots<v_{n}=1\\}$, tolerance $\epsilon$, constraint term
$\mathbf{C}$

2:a set of control points $\\{\mathbf{P}_{kl}\\}_{k=0,l=0}^{n_{1},n_{2}}$
interpolating the point set $\\{\mathbf{Q}_{ij}\\}^{m_{1},m_{2}}_{i=0,j=0}$

3:$\\{\mathbf{P}_{kl}\\}_{k=0,l=0}^{n_{1},n_{2}}\leftarrow\\{\mathbf{Q}_{ij}\\}^{m_{1},m_{2}}_{i=0,j=0}$
$\triangleright$ Initialize the control points with sample points

4:$FP\leftarrow$EvaluateFittedSurfacePoints($\\{\mathbf{P}_{kl}\\}_{k=0,l=0}^{n_{1},n_{2}},U,V,$

5:$\\{\xi_{i},\zeta_{j}\\}^{m_{1},m_{2}}_{i=0,j=0},w_{ij}$)

6:$\delta_{ij}\leftarrow$
AddConstraints($\\{\mathbf{P}_{kl}\\}_{k=0,l=0}^{n_{1},n_{2}},\mathbf{C}$)

7:$\Delta_{ij}\leftarrow$
CalculateDeviation($FP,\\{\mathbf{Q}_{ij}\\}^{m_{1},m_{2}}_{i=0,j=0},\delta_{ij}$)

8:while $\Delta_{ij}>\epsilon$ do

9: $\\{\mathbf{P}_{kl}\\}_{k=0,l=0}^{n_{1},n_{2}}\leftarrow$
RegenerateControlPoints($\Delta_{ij}$)

10:
$FP\leftarrow$EvaluateFittedSurfacePoints($\\{\mathbf{P}_{kl}\\}_{k=0,l=0}^{n_{1},n_{2}},U,V,$

11:$\\{\xi_{i},\zeta_{j}\\}^{m_{1},m_{2}}_{i=0,j=0},w_{ij}$)

12: $\delta_{ij}\leftarrow$
AddConstraints($\\{\mathbf{P}_{kl}\\}_{k=0,l=0}^{n_{1},n_{2}},\mathbf{C}$)

13: $\Delta_{ij}\leftarrow$
CalculateDeviation($FP,\\{\mathbf{Q}_{ij}\\}^{m_{1},m_{2}}_{i=0,j=0},\delta_{ij}$)

14:end while

###  3.5 CPIA convergence proof

Considering the page limit, only the convergence proof for Gyroid CPIA is
given here. Those for Diamond and Schwarz_P basically follow the same proof
procedures but with different equations. We have put them in Section 2 of the
supplementary material.

######  Proposition 1.

The CPIA method for Gyroid is convergent and the limit surface is the least-
square approximation outcome of the initial data
$\mathrm{\\{\mathbf{Q}_{ij}\\}^{m_{1},m_{2}}_{i=0,j=0}}$.

###### Proof.

Using CPIA, a sequence of control points and offset surfaces is generated. To
show its convergence, let
$\mathbf{P}^{k}=\\{\mathbf{P}_{0}^{k},\mathbf{P}_{1}^{k},\cdots,\mathbf{P}_{n}^{k}\\}^{T}$
and $\mathbf{Q}=\\{\mathbf{Q}_{0},\mathbf{Q}_{1},\cdots,\mathbf{Q}_{m}\\}^{T}$
denote the control points at $k$-th iteration and the sample points,
respectively. In the $(k+1)$th iteration, we have:

$\displaystyle\begin{aligned}
\mathbf{P}^{k+1}=&\mathbf{P}^{k}+(\mathbf{Q-BwP}^{k})\\\
&+\frac{1}{2}\sum_{i=1}^{4}((\mathbf{M}_{2ig}-\mathbf{M}_{1ig})\mathbf{P}^{k}\mathbf{T}_{(i+3)g}-(\mathbf{M}_{2ig}-\mathbf{M}_{1ig})\mathbf{P}^{k})\\\
&+\frac{1}{2}\sum_{i=1}^{4}((3\mathbf{M}_{4ig}-2\mathbf{M}_{3ig}-\mathbf{M}_{5ig})\mathbf{P}^{k}\mathbf{T}_{(i+3)g}\\\
&-(3\mathbf{M}_{4ig}-2\mathbf{M}_{3ig}-\mathbf{M}_{5ig})\ \mathbf{P}^{k})\\\
\end{aligned}$ (18)

The matrices $\mathbf{T}_{ig}$ $(i=4,5,6,7)$ are mentioned in Sec. 3.4, and
they are invertible matrices where the absolute values of the eigenvalues are
all 1, i.e., $\left|\lambda_{i}(\mathbf{T}_{jg})\right|=1,i=1,2,3,4$ and
$j=4,5,6,7$. The matrices $\mathbf{M}_{ijg}$ $(i=1,2,3,4,5$ and $j=1,2,3,4)$
are also given in Sec. 3.4. $\mathbf{B}$ is the B-spline basis function matrix
and $\mathbf{w}$ is the weight matrix. They are invertible. Let
$\mathbf{r}^{k}=\mathbf{P}^{k}-\mathbf{w}^{-1}\mathbf{B}^{-1}\mathbf{Q}$, $\
\mathbf{D}=\mathbf{I}-\mathbf{Bw}$,
$\alpha_{i}=3\mathbf{M}_{4ig}-2\mathbf{M}_{3ig}-\mathbf{M}_{5ig}$, and
$\beta_{i}=\mathbf{M}_{2ig}-\mathbf{M}_{1ig}(i=1,2,3,4)$. Then Eq. (18) can be
rewritten as:

$\displaystyle\begin{aligned}
&\mathbf{r}^{k+1}=(\mathbf{D}-\frac{1}{2}\sum_{i=0}^{4}(\alpha_{i}+\beta_{i}))\mathbf{r}^{k}+\frac{1}{2}\sum_{i=0}^{4}(\alpha_{i}+\beta_{i})\mathbf{r}^{k}\mathbf{T}_{(i+3)g}\\\
&=(\mathbf{D}-\frac{1}{2}\sum_{i=0}^{4}(\alpha_{i}+\beta_{i}))((\mathbf{D}-\frac{1}{2}\sum_{i=0}^{4}(\alpha_{i}+\beta_{i}))\mathbf{r}^{k-1}+\frac{1}{2}\sum_{i=0}^{4}(\alpha_{i}+\beta_{i})\mathbf{r}^{k-1}\mathbf{T}_{(i+3)g})\\\
&+\frac{1}{2}\sum_{i=0}^{4}(\alpha_{i}+\beta_{i})((\mathbf{D}-\frac{1}{2}\sum_{n=0}^{4}(\alpha_{n}+\beta_{n}))\mathbf{r}^{k-1}+\frac{1}{2}\sum_{n=0}^{4}(\alpha_{n}+\beta_{n})\mathbf{r}^{k-1}\mathbf{T}_{(n+3)g})\mathbf{T}_{(i+3)g}\\\
&=\cdots\\\ &=\sum_{i=0}^{k+1}{k+1\choose
i}(\mathbf{D}-\frac{1}{2}\sum_{n=0}^{4}(\alpha_{n}+\beta_{n}))^{i}(\frac{1}{2}\sum_{n=0}^{4}(\alpha_{n}+\beta_{n}))^{k+1-i}\mathbf{r}^{0}\mathbf{T}_{(i+3)g}^{k+1-i}\end{aligned}$
(19)

Supposing $\\{\lambda_{k}(\alpha_{i})\\}(k=0,1,\cdots,n-1)$ and
$\\{\lambda_{l}(\beta_{j})\\}(l=0,1,\cdots,n-1)$ are the eigenvalues of
$\alpha_{i}$ and $\beta_{j}$ sorted in non-decreasing order. n is the number
of control points of the approximation surface in the u and v directions.
Since $\mid\lambda_{k}(\alpha_{i})\mid=0\quad or\quad 1,k=0,1,\cdots,n-1$,
$\mid\lambda_{l}(\beta_{j})\mid=0\quad or\quad 1,l=0,1,\cdots,n-1$, and
$r(\sum_{i=0}^{4}\beta_{i})<n$, their powers have no effect on convergence.
$r(\sum_{i=0}^{4}\beta_{i})$ is the rank of $(\sum_{i=0}^{4}\beta_{i})$. Let

$\displaystyle\mathbf{s}^{k+1}=\sum_{i=0}^{k+1}{k+1\choose
i}(\mathbf{D}-\frac{1}{2}\sum_{n=0}^{4}(\alpha_{n}+\beta_{n}))^{i}(\frac{1}{2}\sum_{n=0}^{4}(\alpha_{n}+\beta_{n}))^{k+1-i}\mathbf{r}^{0}$
(20)

To show that $\\{\mathbf{r}^{k}\\}$ is convergent, we firstly show that
$\\{\mathbf{s}^{k}\\}$ is convergent. With the Binomial theorem, we get:

$\mathbf{s}^{k+1}=\mathbf{D}^{k+1}\mathbf{r}^{0}=(\mathbf{I}-\mathbf{Bw})^{k+1}\mathbf{r}^{0}$
(21)

From Theorem 2.2 in [41], we know that $\rho(\mathbf{I}-\mathbf{B})<1$, where
$\rho(\cdot)$ gets a matrix’s spectral radius. With uniform weight assignment,
$\rho(\mathbf{I}-\mathbf{Bw})<1$. Therefore, we have the following equation:

$\displaystyle\lim_{k\to\infty}(\mathbf{I}-\mathbf{Bw})^{k}=\lim_{k\to\infty}\mathbf{s}^{k}=(\mathbf{0})_{n+1}$
(22)

Let

$\displaystyle\begin{aligned}
\mathbf{r}^{k+1}&=\sum_{i=0}^{k+1}\mathbf{r}_{i}=\sum_{i=0}^{k+1}{k+1\choose
i}(\mathbf{D}-\frac{1}{2}\sum_{n=0}^{4}(\alpha_{n}+\beta_{n}))^{i}(\frac{1}{2}\sum_{n=0}^{4}(\alpha_{n}+\beta_{n}))^{k+1-i}\mathbf{r}^{0}\mathbf{T}_{(i+3)g}^{k+1-i}\\\
\mathbf{s}^{k+1}&=\sum_{i=0}^{k+1}\mathbf{s}_{i}=\sum_{i=0}^{k+1}{k+1\choose
i}(\mathbf{D}-\frac{1}{2}\sum_{n=0}^{4}(\alpha_{n}+\beta_{n}))^{i}(\frac{1}{2}\sum_{n=0}^{4}(\alpha_{n}+\beta_{n}))^{k+1-i}\mathbf{r}^{0}\\\
\end{aligned}$ (23)

For each term $\mathbf{r}_{i}$ in $\mathbf{r}^{k+1}$ and $\mathbf{s}_{i}$ in
$s^{k+1}$, $||\mathbf{r}_{i}||_{2}=||\mathbf{s}_{i}||_{2}$. $||\cdot||_{2}$
gets 2-norm of a matrix. Then
$||\mathbf{r}^{k+1}||_{2}=||\mathbf{s}^{k+1}||_{2}$ since
$\mathbf{T}_{ig}(i=4,5,6,7)$ are rigid transformation matrices with modulo
lengths $1$ for all the eigenvalues. Because
$||\mathbf{r}^{k+1}||_{2}=\sqrt[]{\lambda_{max}((\mathbf{r}^{k+1})^{T}\mathbf{r}^{k+1})}$
and
$||\mathbf{s}^{k+1}||_{2}=\sqrt[]{\lambda_{max}((\mathbf{s}^{k+1})^{T}\mathbf{s}^{k+1})}$,
the spectral radius of $\mathbf{r}^{k+1}$ and $\mathbf{s}^{k+1}$ are the same.
So if $\\{\mathbf{s}^{k}\\}$ is convergent, then $\\{\mathbf{r}^{k}\\}$ is
also convergent, which is equivalent to:

$\displaystyle\lim_{k\to\infty}\mathbf{P}^{k}-\mathbf{w}^{-1}\mathbf{B}^{-1}\mathbf{Q}=\mathbf{0}$
(24)

So { $\mathbf{P}^{k}$ } is convergent, and

$\displaystyle\mathbf{P}^{\infty}=\mathbf{w}^{-1}\mathbf{B}^{-1}\mathbf{Q}$
(25)

which concludes the proof. ∎

###  3.6 STEP file generation

As is shown in Fig. 6, a STEP file is divided into two sections. The first is
the head and tail of the file, which contains the meta-data like the author,
the file name, and the description of the file. The second contains all the
data of the file, including geometry information (e.g., “CARTESIAN_POINT”),
topological information (e.g., “CLOSED_SHELL”), and other parameters such as
the tolerance of the model (e.g., “UNCERTAINTY_MEASURE_WITH_UNIT”). In this
region, Express language is used. Express is a formal language designed for
storing and transmitting information about engineering data and processes in
computer systems. In this paper, the TPMS solid model data is stored in STEP
files with this language.

Figure 6: An example of the structure of a STEP file. The structure of a STEP
file includes two regions: the head and tail region and the data region. The
geometry data, topological data, and other data in data regions are framed in
the graph. Figure 7: Illustration of gap patches and offset patches: (a)
Copied and transformed offset patches without generating gap patches; and (b)
Offset patches and corresponding gap patches.

The STEP file generation process used in this work is a standard one and
consists of three stages. First, meta-data is written. Second, the NURBS
patches are copied, transformed, and assembled into a solid model. Third, the
data of the solid model is translated into the Express language and saved in
STEP files.

The above three stages are quite straightforward, except for one special task
before assembly in the second stage. That is, gap patches need to be crafted
to bridge the space between a pair of offset surfaces, as shown in Fig. 7. To
fill the gaps, we first identify pairs of opposite edges on offset surfaces
and then generate surfaces connecting those edge pairs by using the bridging
operation offered by existing CAD software.

Figure 8: Gyroid (left), Diamond (middle), and Schwarz_P (right) models and
analysis results for case study (1)-(6). (a)-(c) unit models and assembled
models; (d)-(f) variation of the second-order partial derivative; (g)-(i)
scaled models. Figure 9: Solid models (left), transparent models (middle), and
the section of the models (right) for case study (7)-(9): the top row for
models designed with Gyroid; the middle row for models designed with Diamond;
and the bottom row for models designed with Schwarz_P. Figure 10: Error
distribution in the NURBS surfaces with five different offset values for error
analysis: the top row for Gyroid; the middle row for Diamond; and the bottom
for Schwarz_P.

With gap patches in place, they are copied and arranged alongside
corresponding offset patches according to the intrinsic symmetry of TPMS. In
this process, if we first create a unit (with gaps filled already) and then
duplicate it alongside the given directions, overlapped gap patches will cause
incorrect model topology. To avoid this problem, we opted to not fill the gaps
immediately after creating a unit, but to keep this unit open. Then the open
unit is duplicated alongside the given directions, resulting in two shells.
The gap between these two shells is then filled. This way of working can avoid
overlapping gap patches inside the solid model, and thus there is no need for
deleting gap patches.

Then all the patches (including offset patches and gap patches) are assembled
to form a TPMS unit (shown in Fig. 7). Due to the spatially periodic nature of
TPMS, the assignment of spatial periods enables precise transformation and
assembly of the patches into a complete solid model. Furthermore, because most
surface patches share the same geometry and they only differ in their
positions and orientations, we only need to store geometry once and all
topological data can be concisely stored by referencing this geometry,
accompanied by their rigid transformation matrices. This trick can
significantly reduce the STEP file for a TPMS model, especially for those
having many cells.

##  4 Results and discussions

The proposed method has been implemented using C++ on a computer with an Intel
Core i9-12900K CPU and 128GB RAM. Based on this implementation, nine case
studies are to be presented to demonstrate the effectiveness of the proposed
method. Case studies 1-3 (Fig. 8a-c) considered three simple uniform
arrangements of typical TPMS structures, including Gyroid, Diamond, and
Schwarz_P. Case studies 4-6 (Fig. 8g-i) analyzed a more complex situation
where TPMS structures were scaled to have distorted shapes, which increased
the difficulty of surface approximation. Case studies 7-9 (Fig. 9) added even
more complexity with trimmed TPMS structures. Error analysis (Figs. 10-13 and
Table 2) of the translated models in those case studies are given. Continuity
analysis (Table 3) and efficiency analysis (Figs. 14 and 15) are also
provided. To test the validity of .step files after translation, Siemens NX
(version 2306) has been used to open those files and then carry out model
edits (e.g., boolean operations), and all case studies passed the test.

Figure 11: The statistical results of the error range distribution of Gyroid
for error analysis. Figure 12: The statistical results of the error range
distribution of Diamond for error analysis. Figure 13: The statistical results
of the error range distribution of Schwarz_P for error analysis. Figure 14:
Time used in constructing NURBS surfaces under different tolerances for
efficiency analysis. Figure 15: Error reduction in CPIA algorithm for
efficiency analysis.

###  4.1 Examples

Case studies 1-3 (Fig. 8 a-c) considered a TPMS unit model and an assembled
model of Gyroid, Diamond, and Schwarz_P TPMS respectively. Case studies 4-6
(Fig. 8 g-i) demonstrate the scaled TPMS solid models with distorted shapes.
Case study 7-9 (Fig. 9) demonstrates the trimmed TPMS structures. Figs. 8 and
9 show the resulting surfaces/solids, and Table 4 summarizes the numbers of
control points of those approximate surfaces.

Error analysis (Figs. 10-13 and Table 2) considered the error distribution in
translated models in the case studies. It also gives a comparison between the
maximum error of the solid model and the given tolerance. From the statistics,
the maximum errors are all below the given tolerances, which validates the
error-controlling feature of our method.

Table 2: The offset, tolerance, and corresponding maximum errors between the TPMS solid model and the original TPMS model for error analysis. Offset | Tolerance | Maximum Errors
---|---|---
Gyroid | Diamond | Schwarz_P
0 | 0.005 | 0.00024 | 0.0018 | 0.00022
0.1 | 0.01 | 0.00077 | 0.0016 | 0.00041
0.2 | 0.01 | 0.0039 | 0.0040 | 0.00025
0.3 | 0.01 | 0.0044 | 0.0032 | 0.00032
0.4 | 0.01 | 0.0032 | 0.0034 | 0.00063
Table 3: The relative maximum and average differences of the first-order and second-order derivatives between two connected edges of two adjacent approximation surfaces for continuity analysis. Model Type | Offset | Maximum Differences | Average Differences
---|---|---|---
First Order | Second Order | First Order | Second Order
Gyroid | 0 | 0.08% | 0.25% | 0.04% | 0.18%
Gyroid | 0.3 | 0.12% | 0.60% | 0.06% | 0.40%
Diamond | 0 | 0.06% | 0.15% | 0.09% | 0.13%
Diamond | 0.3 | 0.07% | 0.15% | 0.09% | 0.13%
Schwarz_P | 0 | 0.08% | 0.82% | 0.04% | 0.42%
Schwarz_P | 0.3 | 0.09% | 0.80% | 0.04% | 0.44%
Table 4: The numbers of control points in approximation surfaces of case studies (1)-(9). Model Type | Offset | Number of Sample Points | Number of Control Points
---|---|---|---
Gyroid | 0 | 24025 | 24025
0.1 | 10609 | 10609
0.2 | 10609 | 10609
0.3 | 10609 | 10609
0.4 | 10609 | 10609
Diamond | 0 | 2601 | 2601
0.1 | 2601 | 2601
0.2 | 2601 | 2601
0.3 | 2601 | 2601
0.4 | 2601 | 2601
Schwarz_P | 0 | 2601 | 2601
0.1 | 2601 | 2601
0.2 | 2601 | 2601
0.3 | 2601 | 2601
0.4 | 2601 | 2601
|  |  |

Continuity analysis (Table 3) analyzed the continuity between adjacent
approximation surfaces in translated models in the case studies. The relative
maximum and average differences of the first-order and second-order
derivatives between two connected edges of two adjacent surfaces are
calculated to validate the $C^{2}$ continuity-preserving feature.

Efficiency analysis (Figs. 14 and 15) demonstrated the construction time of
error-controlled NURBS surfaces with different tolerances to prove the
efficiency of the proposed method and showed the error reduction process of
CPIA to verify the convergence of the CPIA.

###  4.2 Discussion and Limitations

In all of the above TPMS solid models, the NURBS surfaces are constructed with
an error less than the specified tolerance, confirming the effectiveness of
the error-driven TPMS sampling algorithm. The high-accuracy NURBS surface with
an error upper bound of 1mm could be constructed within 10.7s. From Fig. 10,
the error is mainly found in specific parts of the surface, which may be
influenced by the sampling strategy. From Fig. 11, Fig. 12, and Fig. 13, the
proportion of these regions with large errors is quite low. In Fig. 15, the
reduced error during CPIA shows the convergence of CPIA from an experimental
perspective and further confirms the effectiveness of the convergence proof in
Sec. 3.5.

From the error stats in Table 2, we can find that the error bound $2\epsilon$
derived in Section 3.2 is somewhat conservative. For example, the empirical
maximum error is 0.0044 in the case with $offset=0.3$ for Gyroid, but the
theoretical tolerance is 0.01. For Schwarz_P with $offset=0.3$, the maximum
error is only $3.2\%$ of the tolerance, which implies that a tighter bound may
be derived in future work.

In Fig. 8d-f, the variation of the second-order partial derivative
demonstrates the continuous distribution of the second-order partial
derivatives, and in Table 3, the relative differences of the first-order and
second-order derivatives between the connected edges of two adjacent surfaces
are no more than 0.12% and 0.82%, respectively. They both show the $C^{2}$
continuity feature of our method.

In this paper, although only the three most commonly used TPMS types are
considered, this method could also be applied to part of other types of TPMS
structures like I-WP or Schwarz_S. They have different Weierstrass functions
from the three types of TPMS structures utilized in this paper [49]. However,
the Weierstrass functions cannot cover all TPMS models (e.g., the
$\mathbf{C}(\mathbf{P})$ surface) [43]. This method cannot be applied to TPMS
structures with unknown Weierstrass functions. This is a serious limitation of
the proposed method.

##  5 Conclusion

A new method has been presented in this paper to translate TPMS to STEP. The
main features of this method include error control and continuity preservation
in the translation. These main features are essentially achieved through two
new algorithms: an approximation error-driven TPMS sampling algorithm and a
CPIA algorithm. When working together, they can give an error-controlled,
$C^{2}$ continuity-preserving, and fast translation of TPMS models to STEP
files. Both the theoretical proof of the proposed method’s convergence and the
empirical validation of the proposed method’s effectiveness have been
demonstrated by using a series of examples and comparisons.

Although the proposed method is seen to be quite effective in the case studies
and analyses conducted, there are still a few limitations. Certain types of
TPMS cannot be represented by Weierstrass functions (e.g., the
$\mathbf{C}(\mathbf{P})$ surface) [43], thus falling outside the applicability
scope of the proposed method. Another limitation is that the CPIA algorithm
requires the same number of control points as the sample points, leading to a
large number of control points when a high TPMS translation accuracy is
required. Extending the present work to overcome this limitation can be very
practically beneficial, and are among the CAD research studies to be carried
out in the future. The LSPIA algorithm [17] may help here, but further
development is needed to make its approximation error controllable. Surface
approximation with parameter correction, e.g., those presented in [50, 51],
can be helpful in attaining lower approximation errors for PIA and is among
our future research studies.

## Acknowledgments

This work has been funded by NSF of China (No. 62102355), the “Pioneer” and
“Leading Goose” R&D Program of Zhejiang Province (No. 2024C01103), NSF of
Zhejiang Province (No. LQ22F020012), and the Fundamental Research Funds for
the Zhejiang Provincial Universities (No. 2023QZJH32).

## References

## References

  * [1] J. Feng, J. Fu, Z. Lin, C. Shang, X. Niu, Layered infill area generation from triply periodic minimal surfaces for additive manufacturing, Computer-Aided Design 107 (2019) 50–63.
  * [2] S. Wang, Y. Jiang, J. Hu, X. Fan, Z. Luo, Y. Liu, L. Liu, Efficient representation and optimization of TPMS-based porous structures for 3D heat dissipation, Computer-Aided Design 142 (2022) 103123.
  * [3] J. Ding, Q. Zou, S. Qu, P. Bartolo, X. Song, C. C. Wang, STL-free design and manufacturing paradigm for high-precision powder bed fusion, CIRP Annals 70 (2021) 167–170.
  * [4] F. Careri, R. H. Khan, C. Todd, M. M. Attallah, Additive manufacturing of heat exchangers in aerospace applications: a review, Applied Thermal Engineering 235 (2023) 121387.
  * [5] N. Yang, Z. Quan, D. Zhang, Y. Tian, Multi-morphology transition hybridization CAD design of minimal surface porous structures for use in tissue engineering, Computer-Aided Design 56 (2014) 11–21.
  * [6] S. Catchpole-Smith, R. Sélo, A. Davis, I. Ashcroft, C. Tuck, A. Clare, Thermal conductivity of TPMS lattice structures manufactured via laser powder bed fusion, Additive Manufacturing 30 (2019) 100846.
  * [7] Q. Y. Hong, G. Elber, M.-S. Kim, Implicit functionally graded conforming microstructures, Computer-Aided Design 162 (2023) 103548.
  * [8] Q. Zou, H.-Y. Feng, S. Gao, Variational direct modeling: A framework towards integration of parametric modeling and direct modeling in CAD, Computer-Aided Design 157 (2023) 103465.
  * [9] P. Liu, B. Sun, J. Liu, L. Lu, Parametric shell lattice with tailored mechanical properties, Additive Manufacturing 60 (2022) 103258.
  * [10] Y. Xu, H. Pan, R. Wang, Q. Du, L. Lu, New families of triply periodic minimal surface-like shell lattices, Additive Manufacturing 77 (2023) 103779.
  * [11] D.-J. Yoo, Computer-aided porous scaffold design for tissue engineering using triply periodic minimal surfaces, International Journal of Precision Engineering and Manufacturing 12 (1) (2011) 61–71.
  * [12] J. Iamsamang, P. Naiyanetr, Computational method and program for generating a porous scaffold based on implicit surfaces, Computer Methods and Programs in Biomedicine 205 (2021) 106088.
  * [13] M. S. Flores-Jimenez, A. Delgado-Gutiérrez, R. Q. Fuentes-Aguilar, D. Cardenas, Generation of a quadrilateral mesh based on NURBS for gyroids of variable thickness and porosity, Journal of Applied and Computational Mechanics 8 (2) (2022) 684–698.
  * [14] Q. Zou, H.-Y. Feng, Push-pull direct modeling of solid CAD models, Advances in Engineering Software 127 (2019) 59–69.
  * [15] W. Xiao, Y. Liu, R. Li, W. Wang, J. Zheng, G. Zhao, Reconsideration of T-spline data models and their exchanges using STEP, Computer-Aided Design 79 (2016) 36–47.
  * [16] B. Wu, J. Ma, L. Wei, X. Liao, J. Lu, NURBS interpolator with scheduling scheme combining cubic and quartic S-shaped feedrate profiles under drive and chord error constraints, Computer-Aided Design 152 (2022) 103380.
  * [17] H. Lin, T. Maekawa, C. Deng, Survey on geometric iterative methods and their applications, Computer-Aided Design 95 (2018) 40–51.
  * [18] C. Hu, H. Lin, Heterogeneous porous scaffold generation using trivariate B-spline solids and triply periodic minimal surfaces, Graphical Models 115 (2021) 101105.
  * [19] L. Makatura, B. Wang, Y.-L. Chen, B. Deng, C. Wojtan, B. Bickel, W. Matusik, Procedural metamaterials: a unified procedural graph for metamaterial design, ACM Transactions on Graphics 42 (5).
  * [20] Y. Jiang, J. Hu, S. Wang, N. Lei, Z. Luo, L. Liu, Meshless optimization of triply periodic minimal surface based two-fluid heat exchanger, Computer-Aided Design 162 (2023) 103554.
  * [21] J. Feng, J. Fu, C. Shang, Z. Lin, B. Li, Sandwich panel design and performance optimization based on triply periodic minimal surfaces, Computer-Aided Design 115 (2019) 307–322.
  * [22] R. Asbai-Ghoudan, S. Ruiz de Galarreta, N. Rodriguez-Florez, Analytical model for the prediction of permeability of triply periodic minimal surfaces, Journal of the Mechanical Behavior of Biomedical Materials 124 (2021) 104804.
  * [23] K. Karčiauskas, J. Peters, Bi-cubic scaffold surfaces, Computer-Aided Design 150 (2022) 103310.
  * [24] N. Qiu, Y. Wan, Y. Shen, J. Fang, Experimental and numerical studies on mechanical properties of TPMS structures, International Journal of Mechanical Sciences 261 (2024) 108657.
  * [25] G. Savio, R. Meneghello, G. Concheri, Design of variable thickness triply periodic surfaces for additive manufacturing, Progress in Additive Manufacturing 4 (3) (2019) 281–290.
  * [26] S. Rosso, A. Curtarello, F. Basana, L. Grigolato, R. Meneghello, G. Concheri, G. Savio, Modeling symmetric minimal surfaces by mesh subdivision, in: Advances on Mechanics, Design Engineering and Manufacturing III, 2021, pp. 249–254.
  * [27] Q. Pan, G. Xu, Construction of minimal subdivision surface with a given boundary, Computer-Aided Design 43 (4) (2011) 374–380.
  * [28] R. Aubry, S. Dey, E. Mestreau, B. Karamete, D. Gayman, A robust conforming NURBS tessellation for industrial applications based on a mesh generation approach, Computer-Aided Design 63 (2015) 26–38.
  * [29] R. Noruzi, S. Ghadai, O. R. Bingol, A. Krishnamurthy, B. Ganapathysubramanian, NURBS-based microstructure design for organic photovoltaics, Computer-Aided Design 118 (2020) 102771.
  * [30] L. Piegl, W. Tiller, The NURBS Book, Springer-Verlag, Berlin, Heidelberg, 1995.
  * [31] J.-J. Fang, C.-L. Hung, An improved parameterization method for B-spline curve and surface interpolation, Computer-Aided Design 45 (6) (2013) 1005–1028.
  * [32] E. Lee, Choosing nodes in parametric curve interpolation, Computer-Aided Design 21 (6) (1989) 363–370.
  * [33] C.-G. Lim, A universal parametrization in B-spline curve and surface interpolation, Computer Aided Geometric Design 16 (5) (1999) 407–422.
  * [34] T. A. Foley, G. M. Nielson, Knot selection for parametric spline interpolation, in: Mathematical Methods in Computer Aided Geometric Design, 1989, pp. 261–271.
  * [35] A. Iglesias, A. Gálvez, M. Collantes, Four adaptive memetic bat algorithm schemes for Bézier curve parameterization, Transactions on Computational Science XXVIII: Special Issue on Cyberworlds and Cybersecurity 9590 (2016) 127–145.
  * [36] L. Piegl, Surface approximation to scanned data, The Visual Computer 16 (2000) 386–395.
  * [37] J. Luo, H. Kang, Z. Yang, Knot calculation for spline fitting based on the unimodality property, Computer Aided Geometric Design 73 (2019) 54–69.
  * [38] D. Michel, A. Zidna, A new deterministic heuristic knots placement for B-Spline approximation, Mathematics and Computers in Simulation 186 (2021) 91–102.
  * [39] G. Farin, From conics to NURBS: A tutorial and survey, IEEE Computer Graphics and Applications 12 (5) (1992) 78–86.
  * [40] N. Carlson, M. Gulliksson, Surface fitting with nurbs: a gauss newton with trust region approach, in: Proceedings of the 13th WSEAS international conference on Applied mathematics, 2008, pp. 169–174.
  * [41] H.-W. Lin, H.-J. Bao, G.-J. Wang, Totally positive bases and progressive iteration approximation, Computers & Mathematics with Applications 50 (3) (2005) 575–586.
  * [42] P. J. Gandy, D. Cvijović, A. L. Mackay, J. Klinowski, Exact computation of the triply periodic D (’diamond’) minimal surface, Chemical Physics Letters 314 (5-6) (1999) 543–551.
  * [43] P. J. Gandy, J. Klinowski, Exact computation of the triply periodic G (’Gyroid’) minimal surface, Chemical Physics Letters 321 (2000) 363–371.
  * [44] P. Gandy, J. Klinowski, Exact computation of the triply periodic Schwarz P minimal surface, Chemical Physics Letters 322 (2000) 579–586.
  * [45] E. Carberry, Chapter 18: Weierstrass-Enneper representations, https://ocw.mit.edu/courses/18-994-seminar-in-geometry-fall-2004/57baa4cfbbefdf7a1ac1776205c3ad1b_chapter18.pdf, access: 2024-04-10.
  * [46] D. Filip, R. Magedson, R. Markot, Surface algorithms using bounds on derivatives, Computer Aided Geometric Design 3 (4) (1987) 295–311.
  * [47] J. Zheng, T. W. Sederberg, Estimating tessellation parameter intervals for rational curves and surfaces, ACM Transactions on Graphics 19 (1) (2000) 56–77.
  * [48] G. Hu, H. Cao, X. Wang, X. Qin, G2 continuity conditions for generalized Bézier-like surfaces with multiple shape parameters, Journal of Inequalities and Applications 2017 (2017) 248.
  * [49] A. Fogden, S. Hyde, Parametrization of triply periodic minimal surfaces. ii. Regular class solutions, Acta Crystallographica Section A 48 (1992) 575–591.
  * [50] Y. Kineri, M. Wang, H. Lin, T. Maekawa, B-spline surface fitting by iterative geometric interpolation/approximation algorithms, Computer-Aided Design 44 (7) (2012) 697–708.
  * [51] P. Bo, X. Mai, W. Meng, C. Zhang, Improving geometric iterative approximation methods using local approximations, Computers & Graphics 116 (C) (2024) 33–45.