# A Low-rank Spline Approximation of Planar Domains

Construction of spline surfaces from given boundary curves is one of the classical problems in computer aided geometric design, which regains much attention in isogeometric analysis in recent years and is called domain parameterization. However, for most of the state-of-the-art parameterization methods, the rank of the spline parameterization is usually large, which results in higher computational cost in solving numerical PDEs. In this paper, we propose a low-rank representation for the spline parameterization of planar domains using low-rank tensor approximation technique, and apply quasi-conformal map as the framework of the spline parameterization. Under given correspondence of boundary curves, a quasi-conformal map with low rank and low distortion between a unit square and the computational domain can be modeled as a non-linear optimization problem. We propose an efficient algorithm to compute the quasi-conformal map by solving two convex optimization problems alternatively. Experimental results show that our approach can produce a bijective and low-rank parametric spline representation of planar domains, which results in better performance than previous approaches in solving numerical PDEs.

Comments

There are no comments yet.

## Authors

• 2 publications
• 2 publications
• ### Volumetric Spline Parameterization for Isogeometric Analysis

Given the spline representation of the boundary of a three dimensional d...
02/02/2019 ∙ by Maodong Pan, et al. ∙ 0

read it

• ### Low-rank signal subspace: parameterization, projection and signal estimation

The paper contains several theoretical results related to the weighted n...
01/24/2021 ∙ by Nikita Zvonarev, et al. ∙ 0

read it

• ### Constructing IGA-suitable planar parameterization from complex CAD boundary by domain partition and global/local optimization

In this paper, we propose a general framework for constructing IGA-suita...
07/03/2017 ∙ by Gang Xu, et al. ∙ 0

read it

• ### Goal-Oriented Adaptive THB-Spline Schemes for PDE-Based Planar Parameterization

This paper presents a PDE-based planar parameterization framework with s...
01/24/2020 ∙ by Jochen Hinz, et al. ∙ 0

read it

• ### Spectral Tensor Train Parameterization of Deep Learning Layers

We study low-rank parameterizations of weight matrices with embedded spe...
03/07/2021 ∙ by Anton Obukhov, et al. ∙ 15

read it

• ### Low Rank Approximation for Smoothing Spline via Eigensystem Truncation

Smoothing splines provide a powerful and flexible means for nonparametri...
11/23/2019 ∙ by Danqing Xu, et al. ∙ 0

read it

• ### Isogeometric computation reuse method for complex objects with topology-consistent volumetric parameterization

Volumetric spline parameterization and computational efficiency are two ...
06/09/2016 ∙ by Gang Xu, et al. ∙ 0

read it

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Given the boundary curves in 2D or 3D, constructing a parametric spline representation to interpolate the given boundary is a fundamental problem in Computer Aided Geometric Design, and Coons surfaces are a classic tool to solve the problem

(farin1999discrete, ). This problem has been revived in recent years due to its applications in isogeometric analysis (IGA), and it is called domain parameterization. Domain parameterization has a great effect on the accuracy and efficiency in subsequent analysis (cohen2010analysis, ; xu2011parameterization, ; pilgerstorfer2014bounding, ). It is a common requirement that the parameterization should be injective, i.e., the mapping from the parametric domain (generally a unit square) to the computational domain is self-intersection free. In addition, the distortion of the map should be as small as possible, i.e., the areas and angles after mapping should be preserved as much as possible. So far many approaches have been proposed to solve the parameterization problem, e.g. the discrete Coons interpolation (farin1999discrete, ), the harmonic mapping (martin2009volumetric, ; nguyen2010parameterization, ; xu2011variational, ), the spring model (gravesen2012planar, ), the nonlinear optimization method (xu2011parameterization, ), parameterization with non-standard B-splines (e.g., T-splines (zhang2012solid, ; zhang2013conformal, ), THB-splines (falini2015planar, )), the method based on Teichmüller mapping (nian2016planar, ) and so on. While all these methods focus on low distortion and bijectivity of the parameterization, the problem of low-rank parameterization is not discussed. In fact, the rank of the parametric spline representation by these methods is usually large, which results in higher computational cost in subsequent isogeometric analysis (mantzaflaris2014matrix, ; juttler2017low, ). Recently, Juetter and his collaborators have observed that reducing the rank of a parameterization can lead to substantial improvements of the overall efficiency of the numerical simulation (mantzaflaris2014matrix, ; mantzaflaris2017low, ; juttler2017low, ). This observation motivates us to explore parameterization techniques which are able to generate low-rank spline representations.

In this paper, by using low-rank tensor approximation technique, we propose a low-rank representation for planar domain parameterization based on quasi-conformal mapping. Quasi-conformal mapping is a natural extension of conformal mapping which preserves angles (lehto1973quasiconformal, ; lui2010compression, ). The angular distortion and the bijectivitiy of a quasi-conformal map can be characterized by a complex function called the Beltrami coefficient (lui2012optimization, ; lui2013texture, ). By optimizing the norm of the Beltrami coefficient and the rank of the spline representation, we are able to find a planar domain parameterization with low rank and low distortion as much as possible.

The remainder of this paper is organized as follows. Section 2 reviews some related work about domain parameterization and the applications of low-rank tensors in science and engineering. Section 3 presents some preliminary knowledge about quasi-conformal mapping and low-rank tensor approximation. In Section 4, we propose a mathematical model followed by an algorithm to compute a low-rank quasi-conformal map for domain parameterization. Section 5 demonstrates some experimental results of our algorithm and its applications in solving numerical PDEs. Comparisons with the state-of-the-art methods are also provided. Finally, we conclude the paper with a summary and future work in Section 6.

## 2 Related work

### 2.1 Domain parameterization

Domain parameterization is one essential step in isogeometric analysis (hughes2005isogeometric, ). The quality of the parameterization greatly influences the numerical accuracy and efficiency of the numerical simulations (cohen2010analysis, ; xu2013optimal, ; pilgerstorfer2014bounding, ). Over the past decade, many approaches have been proposed to solve the parameterization problem. A simple way for domain parameterization is based on discrete Coons patches proposed by Farin and Hansford (farin1999discrete, ). A spring model was suggested by Gravesen et al. to solve the problem  (gravesen2012planar, ). The harmonic functions have many good properties and they were used in (martin2009volumetric, ; nguyen2010parameterization, ; xu2011variational, ) to construct domain parameterizations. These methods are generally computational inexpensive but the resulting parameterization may not be injective—a deficiency that should be avoided in such type of applications. Xu et al. (xu2011parameterization, ) presented a sophisticated nonlinear optimization technique with the injectivity and the quality of the parameterization as an objective. In (falini2015planar, ), THB-splines is used for planar domain parameterization with varying levels of computational complexity. Recently, Nian et al. (nian2016planar, ) proposed an approach for planar domain parameterization based on Teichmüller mapping, which guarantees a bijective and high-quality parameterization. For 3D domains, a framework was developed in (martin2009volumetric, ) to model a single trivariate B-spline from input boundary triangle meshes with genus-zero topology. Aigner et al. (aigner2009swept, ) presented a variational framework for generating NURBS parametrizations of swept volumes, in which the control points can be obtained by solving an optimization problem. Escobar et al. (escobar2011new, ) proposed a solid T-spline modeling algorithm from a surface triangular mesh. Zhang et al. (zhang2012solid, ) developed a mapping-based method to construct rational trivariate solid T-splines for genus-zero geometry from the boundary triangulations. For meshes with more general geometry, they (wang2013trivariate, ) further used the mapping, subdivision and pillowing techniques to generate high quality T-spline representations. In (xu2013constructing, ), the authors proposed a variational harmonic method to construct analysis-suitable parameterization of a computational domain from given CAD boundary information. For models topologically equivalent to a set of cubes and bounded by B-spline surfaces, they (xu2013analysis, ) further studied the volume parameterization of the multi-block computational domain using the nonlinear optimization method proposed in (xu2011parameterization, ). When dealing with more complex geometric shapes, however, single-patch representations do not provide sufficient flexibility. Multi-patch structures are generally constructed to fulfill the task of low distortion parameterization (xu2015two, ; buchegger2017planar, ). In this paper, we focus on 2D domain parameterization.

### 2.2 Applications of low-rank tensor approximation

Low-rank approximation is very helpful for dimension reduction and data compression, and has been successfully applied in many fields like signal processing, computer vision, patter recognition, computer graphics, etc. A thorough survey on this topic is out of the scope of this paper, and we refer the reader to

(markovsky2011low, ; ma2012sparse, ) and references therein. Tensors, as a generation of matrices in higher dimensions have important applications in science and engineering, e.g., psychometrics, psychometrics and data mining (kolda2009tensor, ). The details of low-rank tensor approximation and its applications have been discussed in depth in (grasedyck2013literature, ). Recently, the low-rank tensor optimization has been applied in graphics and geometric modeling community, e.g., in finding the upright orientation of 3D shapes (wang2014upright, ) and in compact implicit surface reconstructions (pan2016compact, ). For other applications of low-rank tensors in geometric modeling and processing, please refer to (xu2015survey, ) and references therein. Juetter and his collaborators recently addressed the problem of low-rank approximation for isogeometric analysis applications. Mantzaflaris et al. (mantzaflaris2014matrix, ) applied low-rank matrix approximation for accelerating the assembly process of stiffness matrices in isogeometric analysis. They further extended their work to 3D case and employed the tensor decomposition technique for Galerkin-based isogeometric analysis, which can reduce the computation time and storage requirements dramatically (mantzaflaris2017low, ). A construction for low-rank tensor-product spline surfaces from given boundary curves is also proposed by Jüttler et al. (juttler2017low, ).

## 3 Preliminaries

In this section, we give some preliminary knowledge about quasi-conformal mapping and low-rank tensor approximation followed by the definition of rank- spline functions.

### 3.1 Quasi-conformal mapping

The most convenient way to explain quasi-conformal mapping is in complex setting. Let be a complex variable with and being the real and imaginary part of respectively, and be the conjugate of , here . For a differentiable complex function , its complex derivatives are defined as and . A complex function defines a map from a complex plane to a complex plane. When , defines a conformal map which preserves angles and maps an infinitesimal circle to an infinitesimal circle. A quasi-conformal map is a generalization of a conformal map which maps an infinitesimal circle to an infinitesimal ellipse.

###### Definition 1.

Suppose is a complex function, where and are two domains in . If is assumed to have continuous partial derivatives, then is quasi-conformal provided it satisfies the Beltrami equation

 ∂f∂¯z=μ(z)∂f∂z (1)

for some complex valued Lebesgue measurable satisfying . is called the Beltrami coefficient of the map .

The Beltrami coefficient determines the angular deviation from conformality. When , the quasi-conformal map becomes conformal. Define the dilatation of at the point by

 K(z)=1+|μ(z)|1−|μ(z)|.

Then a quasi-conformal map takes infinitesimal circles to infinitesimal ellipses with bounded eccentricity given by the dilatation and the orientation of axis rotates an , as shown in Fig. 1. Furthermore, is orientation preserving and bijective provided and .

Besides angular deviation, another quantity that characterizes the area distortion of a map is the Jacobian of the map. In order to eliminate the area difference between the parametric domain (a unit square) and the computational domain, usually scaled Jacobin is employed:

 Js(f)=J(f)AΩ

where is the area of the computational domain .

### 3.2 Low-rank tensor approximation

A tensor is a multidimensional array. More formally, an th-order or -way tensor is an element of the tensor product of vector spaces, each of which has its own coordinate system (kolda2009tensor, ). An th-order tensor is usually denoted by boldface Euler script letters, e.g., . A first-order tensor is a vector, a second-order tensor is a matrix, and tensors of order three or higher are called higher-order tensors.

An th-order tensor is rank one if it can be written as the outer product of vectors, i.e.,

 X=a(1)∘a(2)∘⋯∘a(n),

where denotes the outer product, and .

CP decomposition Let be an th-order tensor, the CP decomposition factorizes into a sum of component rank-one tensors as follows:

 X=R∑r=1λra(1)r∘a(2)r∘⋯∘a(n)r, (2)

where is a positive integer and for . It’s often useful to assume that are normalized to length one with the weights absorbed into .

The rank of a tensor , denoted as , is the smallest number of components in the above expression (2

). The CP decomposition can be considered to be a higher-order generalization of the matrix singular value decomposition (SVD) which can be described as follows:

Singular value decomposition Let be a matrix, the SVD factorizes into a sum of component rank-one matrices as follows:

 X=UΣVH=R∑r=1λrar∘br (3)

where is an

complex unitary matrix,

is an rectangular diagonal matrix with non-negative real numbers on the diagonal, and is an complex unitary matrix. The diagonal entries of are known as the singular values of and are the nonzero singular values. , are the left column vectors of and respectively.

###### Theorem 1.

((eckart1936approximation, )) The best rank- approximation of is given by a truncated SVD of , that is

 ^X=U^ΣVH=k∑r=1λrar∘br (4)

where has a specific rank , and is the same matrix as except that it contains only the largest singular values (the other singular values are replaced by zero). is called rank- approximation of .

From Theorem 1, we can see that the rank of a matrix , denoted as , is equal to the number of nonzero singular values in SVD of . However, is a nonconvex function, and solving a rank-constrained problem is generally NP-hard. Recently several works (recht2010guaranteed, ; cai2010singular, ; liu2013tensor, ) use the trace norm of a matrix to approximately calculate the rank, which leads to a convex optimization problem. The trace norm of is defined as follows

 ∥X∥∗:=∑iσi(X) (5)

where is the th largest singular value of .

### 3.3 Rank-R spline functions

A multivariate function is said to have rank if it can be represented as a sum of separable functions

 f(x1,…,xn)=R∑r=1n∏k=1f(k)r(xk) (6)

where are univariate functions.

Let be an -variate tensor product spline function of -degree () defined over an -dimensional domain :

 g(x1,…,xn)=∑i1,i2,…,inci1,i2,…,inn∏k=1β(k)ik(xk) (7)

where are the control coefficients, are B-spline basis functions in possibly different univariate spline spaces , and each space is defined by a knot vector and a degree . We collect the basis functions in the knot vector

 β(k)(xk)=[β(k)ik(xk)],k=1,…,n,

Let be the -order coefficient tensor associated with the coefficients defined in (7). If the rank of is R and perform the CP decomposition of as

 C(???)=R∑r=1c(1)r∘c(2)r∘⋯∘c(n)r, (8)

then can be expressed in a sum of products

 g(x1,…,xn)=R∑r=1n∏k=1g(k)r(xk) (9)

of the univariate spline functions

 g(k)r(xk)=c(k)r⋅β(k)(xk).

Thus also has rank , and we call is a rank- spline function.

## 4 Parameterization of computational domains via low-rank tensor approximation

### 4.1 Representation of parameterization

Suppose we are given the B-spline representations of the four boundary curves of a computational domain . Our aim is to compute a B-spline representation for the parameterization domain , that is, a map from the unit square to which is bijective, low distortion and low rank. An example is illustrated in Fig. 2.

Assume the parameterization of the computation domain is expressed by a tensor product B-spline function

 P(x,y):=m∑i=0n∑j=0PijMpi(x)Nqj(y) (10)

where are the control points, and are the B-spline basis functions of degree and w.r.t the knot sequences and in respectively. Since we are working in complex settings, we rewrite the parameterization by a complex function

 f(z):=m∑i=0n∑j=0cijMpi(x)Nqj(y), (11)

where , and . Since the boundary curves of the domain are given, is known for , and , .

From the equation (1), the Beltrami coefficient of can be computed as

 μ(f)=(a−b)+√−1(c+d)(a+b)+√−1(c−d) (12)

where

 a=m∑i=0n∑j=0xij∂Mpi(x)∂xNqj(y),b=m∑i=0n∑j=0yijMpi(x)∂Nqj(y)∂y,c=m∑i=0n∑j=0yij∂Mpi(x)∂xNqj(y),d=m∑i=0n∑j=0xijMpi(x)∂Nqj(y)∂y

### 4.2 Low-rank parameterization model

As explained in Section 3.1, the distortion of a quasi-conformal map is determined by its Beltrami coefficient , thus we formulate the parameterization problem as the following model

 argminf∫^Ω|μ(f)|2dz+ω1∫^Ω|∇μ(f)|2dz+ω2rank(C)s.t.∥μ(f)∥∞<1ci0,cin,c0j,cmj (i=0,1,⋯,m,j=0,1,⋯,n) are given (13)

where is a complex matrix whose elements are the coefficients of defined in (11), and are non-negative weights. The first term of the objective function aims to minimize the conformality distortion of , the second term measures the smoothness of and the third term is the low-rank regularization term which tries to reduce the rank of . In terms of the constraints, the first one guarantees that is locally bijective and the second one is the boundary conditions.

### 4.3 Numerical algorithm

Solving the optimization problem (13) for is challenging since it is highly nonlinear and nonconvex. Instead, we set as the auxiliary variable and replace the function with the nuclear norm introduced in Section 3.2. Thus we obtain the following optimization problem

 argminf,ν∫^Ω|ν|2dz+ω1∫^Ω|∇ν|2dz+ω2∥C∥∗s.t.ν=μ(f)∥ν∥∞<1ci0,cin,c0j,cmj (i=0,1,⋯,m,j=0,1,⋯,n) are given (14)

The above problem is relaxed as

 argminf,ν∫^Ω|ν|2dz+ω1∫^Ω|∇ν|2dz+ω2∥C∥∗+ω3∫^Ω|ν−μ(f)|2dzs.t.∥ν∥∞<1ci0,cin,c0j,cmj (i=0,1,⋯,m,j=0,1,⋯,n) are given (15)

For large enough weight , the optimal solution of the model (15) approximates the solution of (14), where is close enough to . To efficiently solve (15), we solve two sub-problems alternatively. More specifically, we set initially. Suppose is obtained at the th iteration. Fixing , we first minimize (15) for to obtain . Then by fixing , we obtain by minimizing (15) for . The procedure runs until for a user-specified . In the following, we will discuss the two sub-problems in detail.

Problem 1 Given , find such that the following objective function is minimized

 argminfω2∥C∥∗+ω3∫^Ω|ν−μ(f)|2dzs.t.ci0,cin,c0j,cmj (i=0,1,⋯,m,j=0,1,⋯,n) are given (16)

Problem (16) is similar to the complex matrix completion problem (cai2010singular, ). However, since is a rational B-spline function, the problem is still hard to solve. Instead we solve the following relaxed model

 argminfω2∥C∥∗+ω3∫^Ω|f¯z−νfz|2dz+λ∥Pr(C)−y∥2 (17)

where is a large positive weight, is the vector whose elements are comprised of , (), and is a linear operator that shapes the boundary elements of into a vector, i.e., ,. Now (17) is a convex optimization problem which can be solved by the alternating direction method of multipliers (ADMM) efficiently. The ADMM can be viewed as an attempt to blend the benefits of dual decomposition and augmented Lagrangian methods, and is used to solve constrained optimization problems with separable objective functions. The basic approach is outlined as follows.

Variable splitting Since the objective function in (17) is the sum of three functions and one of which is dependent on the others, using variable splitting technique leads to the following constrained optimization problem

 argminfω2∥Z∥∗+ω3∫^Ω|f¯z−νfz|2dz+λ∥Pr(C)−y∥2s.t.c=z (18)

where is the vectorization of , i.e., , is an auxiliary matrix of the same size as , and is the vectorization of .

Augmented Lagrangian One typical way for solving (18) is to use the augmented Lagrangian scheme. In our problem, the augmented Lagrangian function is defined as

 Lρ(c,Z,η)=ω2∥Z∥∗+ω3∫^Ω|f¯z−νfz|2dz+λ∥Pr(C)−y∥2+Re(ηH(c−z))+ρ2∥c−z∥2 (19)

where is the real part of the complex number z, is a vector of Lagrangian multiplier corresponding to the constraint , and is the penalty parameter. Now the ADMM algorithm can be outlined as follows

-subproblem  The subproblem for is

 argmincLρ(c,Zt,ηt)=ω3∫^Ω|f¯z−νfz|2dz+λ∥Pr(C)−y∥2+Re((ηt)H(c−zt))+ρ2∥c−zt∥2 (20)

This is a quadratic optimization problem and the solution can be obtained by solving a sparse and symmetric linear system of equations. The preconditioned conjugate gradient method with incomplete Cholesky factorization is applied in our algorithm.

-subproblem  The subproblem for is

 argminZLρ(ct+1,Z,ηt)=ω2∥Z∥∗+Re((ηt)H(ct+1−z))+ρ2∥ct+1−z∥2 (21)

which has the following closed form solution (cai2010singular, ):

 Zt+1=proxtrω2/ρ(ct+1+ηt/ρ), (22)

Note that the argument must be converted into a matrix of the same size as . Here the proximal operator can be considered as a shrinkage operation on the singular values and is defined as follows

 proxtrω2/ρ(Y)=Umax(S−ω2I/ρ,0)VH, (23)

where is the singular value decomposition (SVD) of , and the max operation is taken element-wise. Please refer to (cai2010singular, ) for the detailed derivation.

In our implementation, and are set as zero, and the stopping criterion is that the value of has small change or the maximum number of iterations reaches.

Problem 2 Given a mapping from to , can be computed by (12), and the problem (15) reduces to the following model

 argminν∫^Ω|ν|2dz+ω1∫^Ω|∇ν|2dz+ω3∫^Ω|ν−μ(f)|2dzs.t.∥ν∥∞<1 (24)

Let

 ν=~m∑i=0~n∑j=0~cijMpi(x)Nqj(y), (25)

where , and and are the B-spline basis functions defined in (10). For the simplicity of computation, the constraint in the above optimization problem is replaced by

 −√22<~xij<√22, −√22<~yij<√22,i=0,1,⋯,~m,j=0,1,⋯,~n (26)

Then (24) becomes a quadratic optimization problem which can be easily solved.

### 4.4 Post-processing

The above algorithm iteratively solves two sub-problems to obtain two sequences of complex functions and . In order to accelerate the convergence of the algorithm, we add a weight into the second term of problem (16) after iterations, where satisfies for a threshold , which leads to the following problem

 argminfω2∥C∥∗+ω3∫^Ωω|ν−μ(f)|2dzs.t.ci0,cin,c0j,cmj (i=0,1,⋯,m,j=0,1,⋯,n) are given (27)

where the weight and is a threshold which helps to avoid division by zero. The problem (27) can be solved in the same way as the problem (16). From the numerical examples, we can see that this post-processing step is essential and effective, see Fig. 3 for a comparison result.

Now the overall algorithm of our parameterization method is summarized in Algorithm 2.

###### Remark 1.

Our mathematical model (13) and the registration model presented in (lam2014landmark, ) both obtain a diffeomorphism via quasi-conformal mapping. However, there are several differences between the two methods. Firstly, we not only want to find a quasi-conformal map with low distortion but also add a low-rank regularization term in (13) to make the rank of the map as low as possible. Secondly, in (lam2014landmark, ), the map is represented in a discrete form, while it is expressed in a continuous form, i.e., tensor product B-splines in our work. Finally, the second constraint in the model (13) is ignored in (lam2014landmark, ), which can not guarantee the bijectivity of .

## 5 Results and discussions

In this section, we demonstrate some examples to show the effectiveness of our parameterization method by comparing it with several state-of-the-art methods. The application of our parameterization in solving numerical PDEs is also provided.

### 5.1 Implementation details

We implement our algorithm on a PC with a quad-core Intel i5 @3.1GHz and 8GB of RAM using C++ and MATLAB. There are several parameters for setting. Most of them are set as default values, e.g., the penalty parameter is typically set to be , the threshold and in Section 4.4 are set as and respectively, and the weight in (17) is set to be . We use bicubic uniform B-splines to represent the map and the auxiliary variable (i.e., in (11) and (25)). Unless specified, the knot parameters in (11) and (25) are chosen as and respectively in our examples, which is proven to work well.

There are three weights , and in the mathematical model (15). The weight controls the smoothness of and we typically set . The penalty weight is used to control the difference between and and is set to be 100 in practice. The weight can be used to balance the the rank of the map and parameterization quality. Clearly, larger can reduce the rank of while smaller leads to parameterization results of higher quality. We observe that provides a good compromise between the rank and the quality. Fig. 4 provides an illustrating example.

### 5.2 Parameteriation results

In the following, we present some examples to demonstrate the low rank and low distortion properties of our parameterization method.

#### 5.2.1 Low rank

To demonstrate the superiority of our method in terms of the rank of the map, we provide a comparison with several state-of-the-art parameterization methods: nonlinear optimization method (xu2011parameterization, ), variational harmonic method (xu2013constructing, ), the Teichmüller mapping method (T-Map) (nian2016planar, ) and the low-rank spline interpolation method (low-rank spline) (juttler2017low, ). To have a fair comparison, the number of knots in the B-spline representation (10) are chosen to be the same () for these methods. The rank of the map is numerically computed as the number of singular values of the complex matrix which are greater than a user-specified threshold ( in our experiment). Table 1 shows the statistics of our experiments. Besides the low-rank spline method which sacrifices the parameterization quality, our method significantly outperforms other state-of-the-art parameterization methods in terms of the rank. Some of the parameterization results are shown in Fig 5. As described in Section 3.3, owing to the low-rank property of our method, the map can be represented in a sum of the product of univariate spline functions with a small number of terms, which helps to speed up the assembly process in IGA without sacrificing the overall accuracy of the simulation, see Section 5.2.4 for some examples.

#### 5.2.2 Local injectivity

Fig. 6, Fig. 7 and Fig. 8 depict the parameterization results of the rabbit, the butterfly and the dolphin by different methods respectively. We observe that the variational harmonic method and the low-rank spline method have many self-intersections in the concave regions, the nonlinear optimization method produces non-injective mapping in the butterfly model, the T-map method is not injective in some regions in the Rabbit model, e.g. in the ear of the rabbit, while our method is always injective in these examples.

#### 5.2.3 Distortion

Besides injectivity, the map distortion (including angular distortion described by and the area distortion represented by ) is an important criteria to measure the quality of the parameterization. In our experiments, to measure the area distortion of a map, we firstly uniformly subdivide the parametric domain into sub-rectangles (, ), then the area distortion over the sub-rectangle , denoted as , is calculated as follows

 Js(f)|^Ωij=∫^ΩijJs(f)dxdyA^Ωij (28)

where is the area of .

From Table 1, we can see that the variational harmonic method and the low-rank spline method are significantly worse than our method in terms of distortion. In Fig. 9, Fig. 10, and Fig. 11, we compare our method with the nonlinear optimization method and T-map method by displaying the colormaps of the Beltrami coefficients and the scaled Jacobian . It can be seen that our method produces smaller angular distortion in some concave regions, e.g. in the ear of the rabbit, in the root of the wing of the butterfly, and in the body of the dolphin than the other two methods,  and at the same time, our approach achieves best results in terms of the area distortion among the three methods.

In summary, our method produces much better parameterization results than the other state-of-the-art methods in all these examples. The reason might be as follows. The variational harmonic method and the low-rank spline interpolation method can’t guarantee the injectivity in theory. The Teichmüller mapping method computes a Teichmüller map by solving a nonlinear and non-convex optimization problem. But their method has no convergence guarantee, and thus it may not be able to find the real Teichmüller map in some cases. The nonlinear optimization method puts some strict constraints, which could result in no solutions for complex shapes.

#### 5.2.4 Solving PDEs using IGA

In this subsection, we apply our low-rank parameterization together with IGA to solve numerical partial differential equations (PDEs) on different domains. The stability, accuracy and efficiency of the numerical simulation are compared with the nonlinear optimization method and the T-map method.

Consider the following elliptic problem

 {−Δu+u=finΩu|∂Ω=gon∂Ω (29)

where is a Lipschitz continuous domain with boundary , are given. The variational form of the problem (29) consists in finding , such that

 a(u,v)=f(v),∀v∈H10(Ω). (30)

where

 a(u,v)=∫Ω(∇u⋅∇v+uv)dxdy,f(v)=∫Ωfvdxdy

Let , then the problem (30) is equivalent to find , such that

 a(w,v)=l(v),∀v∈H10(Ω). (31)

where .

In the setting of isogeometric analysis, the domain is parameterized by a global map which is defined in (11). The isogeometric discretization takes advantages of the given parameterization of the domain . In particular, the discretization space can be chosen as

 Vh=span{ϕij(x,y),i=0,1,…,m,j=0,1,…,n}

with and .

The finite-dimensional space is now used for the Galerkin discretization of the variational formulation (31), which consists in finding , such that

 a(wh,vh)=l(vh),∀vh∈Vh. (32)

We solve the above elliptic problem over two domain examples (the Rabbit-shaped domain and the Butterfly-shaped domain) to show the numerical advantages of our low-rank parameterization method.

Rabbit-shaped domain with different parameterizations In this example, we solve the elliptic problem over the Rabbit-shaped domain, where has an exact solution . The parameterization results of this domain are shown in Fig. 6

. The degrees of freedom (

) of the basis functions in in this example is . Fig. 12(a)12(b) and 12(c) show the numerical errors of the solutions for the nonlinear optimization method, the T-map method and our method respectively, and Table 2 summaries the condition numbers of the stiffness matrices, errors and the assembling time for these three methods. We can see that our method produces smaller condition numbers and errors than the other two methods. At the same time, owing to the low-rank property, our method can accelerate the assembly process of the matrices in IGA by using the low-rank assembly strategy presented in (mantzaflaris2014matrix, ; mantzaflaris2016low).

Butterfly-shaped domain with different parameterizations We consider the elliptic problem (29) over another domain—the Butterfly-shaped domain, where the parameterization results of the nonlinear optimization method, T-map method and our method are shown in Fig. 7. The of the basis functions in and the exact solution in this example are and respectively. Fig. 13(a)13(b)13(c) show the numerical errors of the solutions and Table 3 lists the condition numbers of the stiffness matrices, errors and the assembling time for these three methods. Again we can see that our method produces smaller condition numbers and errors than the other two methods, and at the same time, our method can accelerate the assembly process of the stiffness matrices in IGA.

## 6 Conclusions and future work

Parameterization of computational domains and efficiently assembling the mass and stiffness matrices are two essential steps in isogeometric analysis applications. In this paper, using low-rank tensor approximation technique, we propose a low-rank representation scheme for domain parametrization based on quasi-conformal mapping. The problem is formulated as a non-linear and non-convex optimization problem which minimizes the angular distortion and the rank of the map while ensuring the bijectivity of the map. The optimization problem is then converted into two quadratic optimization problems which are solved alternatively. Several experimental examples show that our approach can produce a low-rank and low-distortion parameterization which is superior to other state-of-the-art methods. Numerical examples of our parameterization method together with IGA in solving numerical PDEs also demonstrate some numerical advantages of our method than previous approaches.

Regarding the future work, extending our work to three-dimensional volumetric parametrization is worthy of further research. However, the parametrization problem in three dimensional case is much harder since there is no analogous complex structure in three-dimensional space.

## Acknowledgement

This work is supported by the NSF of China (No. 11571338, 11626253) and by the Fundamental Research Funds for the Central Universities (WK0010000051).

## References

• (1) Aigner M, Heinrich C, J ttler B, et al. Swept Volume Parameterization for Isogeometric Analysis. In: IMA Conference on the Mathematics of Surfaces. 2009: 19-44.
• (2) Buchegger F, Jüttler B. Planar multi-patch domain parameterization via patch adjacency graphs. Computer-Aided Design, 2017, 82: 2-12.
• (3) Cai J F, Cand s E J, Shen Z. A singular value thresholding algorithm for matrix completion. SIAM Journal on Optimization, 2010, 20(4): 1956-1982.
• (4) Cohen E, Martin T, Kirby R M, et al. Analysis-aware modeling: Understanding quality considerations in modeling for isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2010, 199(5): 334-356.
• (5) Eckart C, Young G. The approximation of one matrix by another of lower rank. Psychometrika, 1936, 1(3): 211-218.
• (6) Escobar J M, Casc n J M, Rodr guez E, et al. A new approach to solid modeling with trivariate T-splines based on mesh optimization. Computer Methods in Applied Mechanics and Engineering, 2011, 200(45): 3210-3222.
• (7) Falini A, Špeh J, Jüttler B. Planar domain parameterization with THB-splines. Computer Aided Geometric Design, 2015, 35: 95-108.
• (8) Farin G, Hansford D. Discrete coons patches. Computer Aided Geometric Design, 1999, 16(7): 691-700.
• (9) Grasedyck L, Kressner D, Tobler C. A literature survey of low-rank tensor approximation techniques. GAMM-Mitteilungen, 2013, 36(1): 53-78.
• (10) Gravesen J, Evgrafov A, Nguyen D M, et al. Planar parametrization in isogeometric analysis. In: International Conference on Mathematical Methods for Curves and Surfaces. Springer, Berlin, Heidelberg, 2012: 189-212.
• (11) Hughes T J R, Cottrell J A, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer methods in applied mechanics and engineering, 2005, 194(39): 4135-4195.
• (12) Jüttler B, Mokriš D. Low rank interpolation of boundary spline curves. Computer Aided Geometric Design, 2017.
• (13) Kolda T G, Bader B W. Tensor decompositions and applications. SIAM review, 2009, 51(3): 455-500.
• (14) Lam K C, Lui L M. Landmark-and intensity-based registration with large deformations via quasi-conformal maps. SIAM Journal on Imaging Sciences, 2014, 7(4): 2364-2392.
• (15) Lehto O, Virtanen K I. Quasiconformal mappings in the plane. New York: Springer, 1973.
• (16)

Liu J, Musialski P, Wonka P, et al. Tensor completion for estimating missing values in visual data. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2013, 35(1): 208-220.

• (17) Lui L M, Lam K C, Wong T W, et al. Texture map and video compression using Beltrami representation. SIAM Journal on Imaging Sciences, 2013, 6(4): 1880-1902.
• (18)

Lui L M, Wong T W, Thompson P, et al. Compression of surface registrations using Beltrami coefficients. In: Computer Vision and Pattern Recognition (CVPR), 2010 IEEE Conference on. IEEE, 2010: 2839-2846.

• (19) Lui L M, Wong T W, Zeng W, et al. Optimization of surface registrations using Beltrami holomorphic flow. Journal of scientific computing, 2012, 50(3): 557-585.
• (20) Ma Y, Wright J, Yang A Y. Sparse representation and low-rank representation in computer vision. ECCV Short Course, 2012.
• (21) Mantzaflaris A, Jüttler B, Khoromskij B N, et al. Matrix generation in isogeometric analysis by low rank tensor approximation. In: International Conference on Curves and Surfaces. Springer, Cham, 2014: 321-340.
• (22) Mantzaflaris A, Jüttler B, Khoromskij B N, et al. Low rank tensor methods in Galerkin-based isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2017, 316: 1062-1085.
• (23) Markovsky I. Low rank approximation: algorithms, implementation, applications. Springer Science & Business Media, 2011.
• (24) Martin T, Cohen E, Kirby M. Volumetric parameterization and trivariate B-spline fitting using harmonic functions. In: Proceedings of the 2008 ACM symposium on Solid and physical modeling. ACM, 2008: 269-280.
• (25) Nguyen T, Jüttler B. Parameterization of Contractible Domains Using Sequences of Harmonic Maps. Curves and surfaces, 2010, 6920: 501-514.
• (26) Nian X, Chen F. Planar domain parameterization for isogeometric analysis based on teichm ller mapping. Computer Methods in Applied Mechanics and Engineering, 2016, 311: 41-55.
• (27) Pan M, Tong W, Chen F. Compact implicit surface reconstruction via low-rank tensor approximation. computer-aided design, 2016, 78: 158-167.
• (28) Pilgerstorfer E, Jüttler B. Bounding the influence of domain parameterization and knot spacing on numerical stability in Isogeometric Analysis. Computer Methods in Applied Mechanics and Engineering, 2014, 268: 589-613.
• (29) Recht B, Fazel M, Parrilo P A. Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization. SIAM review, 2010, 52(3): 471-501.
• (30) Wang W, Liu X, Liu L. Upright orientation of 3D shapes via tensor rank minimization[J]. Journal of Mechanical Science and Technology, 2014, 28(7): 2469-2477.
• (31) Wang W, Zhang Y, Liu L, et al. Trivariate solid T-spline construction from boundary triangulations with arbitrary genus topology. Computer-Aided Design, 2013, 45(2): 351-360.
• (32) Xu G, Mourrain B, Duvigneau R, et al. Variational harmonic method for parameterization of computational domain in 2D isogeometric analysis. In: Computer-Aided Design and Computer Graphics (CAD/Graphics), 2011 12th International Conference on. IEEE, 2011: 223-228.
• (33) Xu G, Mourrain B, Duvigneau R, et al. Parameterization of computational domain in isogeometric analysis: methods and comparison. Computer Methods in Applied Mechanics and Engineering, 2011, 200(23): 2021-2031.
• (34) Xu G, Mourrain B, Duvigneau R, et al. Constructing analysis-suitable parameterization of computational domain from CAD boundary by variational harmonic method. Journal of Computational Physics, 2013, 252: 275-289.
• (35) Xu G, Mourrain B, Duvigneau R, et al. Optimal analysis-aware parameterization of computational domain in 3D isogeometric analysis. Computer-Aided Design, 2013, 45(4): 812-821.
• (36) Xu G, Mourrain B, Duvigneau R G, et al. Analysis-suitable volume parameterization of multi-block computational domain in isogeometric applications. Computer-Aided Design, 2013, 45(2): 395-404
• (37) Xu J, Chen F, Deng J. Two-dimensional domain decomposition based on skeleton computation for parameterization and isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 541-555.
• (38) Xu L, Wang R, Zhang J, et al. Survey on sparsity in geometric modeling and processing. Graphical Models, 2015, 82: 160-180.
• (39) Zhang Y, Wang W, Hughes T J R. Solid T-spline construction from boundary representations for genus-zero geometry. Computer Methods in Applied Mechanics and Engineering, 2012, 249: 185-197.
• (40) Zhang Y, Wang W, Hughes T J R. Conformal solid T-spline construction from boundary T-spline representations. Computational Mechanics, 2013: 1-9.