A Novel 3D Mapping Representation and its Applications

In this paper, we focus solely on the left polar decomposition for real rectangular matrices. We begin by denoting the set of $n\times n$ real matrices as $\mathcal{M}_{n}$. The left polar decomposition of a matrix $J\in\mathcal{M}_{n}$ is a factorization of the form $J=UP$, where the columns of the matrix $U\in\mathcal{M}_{n}$ form an orthonormal basis, and $P\in\mathcal{M}_{n}$ is a symmetric positive semi-definite matrix. Specifically, $P=\sqrt{J^{T}J}$.

To better understand this definition, we first note that $J^{T}J$ is a symmetric semi-definite matrix. By applying the eigenvalue decomposition to this matrix, we have

$$J^{T}J=WDW^{-1}$$

where $W$ is an orthogonal matrix because $J^{T}J$ is symmetric, therefore we also have $W^{-1}=W^{T}$. The diagonal matrix $D$ contains the eigenvalues of $J^{T}J$. $\Sigma$ denotes the matrix whose diagonal entries are the square root of the corresponding entries in $D$. Those diagonal entries are the singular values of $J$. Thus, we have $P=\sqrt{J^{T}J}=W\Sigma W^{-1}$.

Therefore, any linear transformation from $\mathbb{R}^{n}$ to $\mathbb{R}^{n}$ is a composition of a rotation/reflection and a dilation.

A surface with a conformal structure is known as a Riemann surface. Given two Riemann surfaces $M$ and $N$, a map $f:M\rightarrow N$ is conformal if it preserves the surface metric up to a multiplicative factor called the conformal factor. Quasiconformal maps are a generalization of conformal maps that are orientation-preserving diffeomorphisms between Riemann surfaces with bounded conformal distortion, meaning that their first-order approximation maps small circles to small ellipses of bounded eccentricity (as shown in Figure 1).

Mathematically, a map $f:\mathbb{C}\rightarrow\mathbb{C}$ is quasiconformal if it satisfies the Beltrami equation

$$\frac{\partial f}{\partial\bar{z}}=\mu(z)\frac{\partial f}{\partial z},$$

for some complex valued functions $\mu$ with $\lVert\mu\rVert_{\infty}<1$ called the Beltrami coefficient. In particular, $f$ is conformal around a small neighborhood of a point $p$ if and only if $\mu(p)=0$, that is

	$\displaystyle f(z)$	$\displaystyle=f(p)+f_{z}(p)z+f_{\bar{z}}(p)\bar{z}$
		$\displaystyle=f(p)+f_{z}(p)(z+\mu(p)\bar{z})$

This equation shows that $f$ can be considered as a map composed of a translation that maps $p$ to $f(p)$, followed by a stretch map $S(z)=z+\mu\bar{z}$, which is postcomposed by a multiplication of a conformal map $f_{z}(p)$. Among these terms, $S(z)$ causes $f$ to map small circles to small ellipses, and $\mu(p)$ controls the stretching effect of $S(z)$. By denoting $z=\rho+i\tau$ and $\mu=re^{i\theta}$, where $r=\left|\mu\right|$, and representing $S(z)=z+\mu\bar{z}$ in a matrix form, we have

$$\left[\begin{pmatrix}1&0\\ 0&1\end{pmatrix}+r\begin{pmatrix}\cos{\theta}&-\sin{\theta}\\ \sin{\theta}&\cos{\theta}\end{pmatrix}\begin{pmatrix}1&0\\ 0&-1\end{pmatrix}\right]\begin{pmatrix}\rho\\ \tau\end{pmatrix}=\begin{pmatrix}1+r\cos{\theta}&r\sin{\theta}\\ r\sin{\theta}&1-r\cos{\theta}\end{pmatrix}\begin{pmatrix}\rho\\ \tau\end{pmatrix}$$

(1)

One can check that the eigenvalues for this matrix are $1+|\mu|$ and $1-|\mu|$. The eigenvectors for this matrix are $(\cos{\frac{\theta}{2}},\sin{\frac{\theta}{2}})^{T}$ and $(\sin{\frac{\theta}{2}},-\cos{\frac{\theta}{2}})^{T}$ respectively. It is worth noting that the two eigenvectors are orthogonal to each other. The distortion or dilation can be measured by $K=\frac{1+|\mu|}{1-|\mu|}$. Thus the Beltrami coefficient $\mu$ gives us important information about the properties of a map.

In addition, we notice some interesting results related to matrix inverse, which facilitate the derivation in the following sections.

Generally, for any invertible matrix $M$, we have

	$\displaystyle M^{-1}$	$\displaystyle=\begin{pmatrix}m_{11}&m_{12}&m_{13}\\ m_{21}&m_{22}&m_{23}\\ m_{31}&m_{32}&m_{33}\\ \end{pmatrix}^{-1}$
		$\displaystyle=\frac{1}{\det(M)}\begin{pmatrix}m_{22}m_{33}-m_{23}m_{32}&m_{13}m_{32}-m_{12}m_{33}&m_{12}m_{23}-m_{13}m_{22}\\ m_{23}m_{31}-m_{21}m_{33}&m_{11}m_{33}-m_{13}m_{31}&m_{21}m_{13}-m_{11}m_{23}\\ m_{21}m_{32}-m_{22}m_{31}&m_{12}m_{31}-m_{11}m_{32}&m_{11}m_{22}-m_{21}m_{12}\end{pmatrix}$

Representing the row vectors of $M$ as $\vec{r}_{1}$, $\vec{r}_{2}$, and $\vec{r}_{3}$, and the column vectors as $\vec{c}_{1}$, $\vec{c}_{2}$, and $\vec{c}_{3}$, we have

	$\displaystyle\begin{pmatrix}\text{---}&\vec{r}_{1}^{T}&\text{---}\\ \text{---}&\vec{r}_{2}^{T}&\text{---}\\ \text{---}&\vec{r}_{3}^{T}&\text{---}\\ \end{pmatrix}^{-1}$	$\displaystyle=\frac{1}{\det(M)}\begin{pmatrix}|&|&|\\ \vec{r}_{2}\times\vec{r}_{3}&\vec{r}_{3}\times\vec{r}_{1}&\vec{r}_{1}\times\vec{r}_{2}\\ |&|&|\end{pmatrix}$		(2)
	$\displaystyle\begin{pmatrix}|&|&|\\ \vec{c}_{1}&\vec{c}_{2}&\vec{c}_{3}\\ |&|&|\end{pmatrix}^{-1}$	$\displaystyle=\frac{1}{\det(M)}\begin{pmatrix}\text{---}&(\vec{c}_{2}\times\vec{c}_{3})^{T}&\text{---}\\ \text{---}&(\vec{c}_{3}\times\vec{c}_{1})^{T}&\text{---}\\ \text{---}&(\vec{c}_{1}\times\vec{c}_{2})^{T}&\text{---}\end{pmatrix}$		(3)

Let $\Omega_{1}$ and $\Omega_{2}$ be oriented, simply connected open subsets of $\mathbb{R}^{3}$. Suppose we have a diffeomorphism $f:\Omega_{1}\rightarrow\Omega_{2}$. $\textbf{J}_{f}(p)$ denotes the Jacobian matrix of $f$ at a point $p$ in $\Omega_{1}$. Since $f$ is diffeomorphic, it is invertible. As a result, the Jacobian determinant $\det(\textbf{J}_{f}(p))$ is always positive for every point $p$ in $\Omega_{1}$.

Let $\Sigma$ be the diagonal matrix containing the singular values $a,b,\text{ and }c$ of $\textbf{J}_{f}(p)$, where by convention, we require that $a\geq b\geq c$. By polar decomposition, we have

$\displaystyle\begin{split}\textbf{J}_{f}(p)&=U(p)P(p)\\ P(p)&=\sqrt{\textbf{J}_{f}(p)^{T}\textbf{J}_{f}(p)}=W\Sigma W^{-1}\end{split}$

(4)

where $U$ is a rotation matrix, because the determinant of $\textbf{J}_{f}(p)$ and $P(p)$ are greater than 0, implying that $\det(U(p))=1$. $a,b\text{ and }c$ are the eigenvalues of $P(p)$ representing the dilation information, and $W$ is an orthogonal matrix containing eigenvectors for dilation.

Recall that in 2D space, the Beltrami coefficient $\mu$ encodes the two dilation $1+|\mu|$ and $1-|\mu|$, together with their corresponding directions $(\cos{\frac{\theta}{2}},\sin{\frac{\theta}{2}})^{T}$ and $(\sin{\frac{\theta}{2}},-\cos{\frac{\theta}{2}})^{T}$, which are orthogonal to each other. Now in 3D space, we have three dilations $a$, $b$, and $c$, together with their corresponding directions represented by the column vectors of $W$, which are also orthogonal to each other, as $P(p)$ is symmetric. Therefore, similar to the Beltrami coefficient, $P(p)$ encodes the dilation information of the 3D mapping $f$.

Furthermore, it is possible to represent the matrix $W$ by Euler angles. To do so, we require that the column vectors of $W$ form an orthonormal basis, and $\det(W)=1$. However, the $\det(W)$ is not necessarily positive, meaning that $W$ may not be a rotation matrix. This can be fixed by randomly multiplying a column of $W$ by $-1$ to get a matrix $\tilde{W}$. This modification ensures $\det(\tilde{W})=1$, indicating that $\tilde{W}$ is a rotation matrix, while at the same time, the matrix $P(p)$ will not be changed. To see this, the matrix $P(p)$ can be reorganized as $a\vec{w}_{1}\vec{w}_{1}^{T}+b\vec{w}_{2}\vec{w}_{2}^{T}+c\vec{w}_{3}\vec{w}_{3}^{T}$, where $\vec{w}_{1},\vec{w}_{2}\text{ and }\vec{w}_{3}$ are the column vectors of $W$. Let $\tilde{\vec{w}}_{1}=-\vec{w}_{1}$, then $\tilde{\vec{w}}_{1}\tilde{\vec{w}}_{1}^{T}=\vec{w}_{1}\vec{w}_{1}^{T}$. Thus reversing the direction of $\vec{w}_{1}$ makes no change to $P(p)$. This also holds for $\vec{w}_{2}$ and $\vec{w}_{3}$.

As is known, a rotation matrix $R$ can be written in the following form:

	$\displaystyle R$	$\displaystyle=R_{z}(\theta_{z})R_{y}(\theta_{y})R_{x}(\theta_{x})$		(5)
		$\displaystyle=\begin{pmatrix}\cos{\theta_{z}}&-\sin{\theta_{z}}&0\\ \sin{\theta_{z}}&\cos{\theta_{z}}&0\\ 0&0&1\end{pmatrix}\begin{pmatrix}\cos{\theta_{y}}&0&\sin{\theta_{y}}\\ 0&1&0\\ -\sin{\theta_{y}}&0&\cos{\theta_{y}}\end{pmatrix}\begin{pmatrix}1&0&0\\ 0&\cos{\theta_{x}}&-\sin{\theta_{x}}\\ 0&\sin{\theta_{x}}&\cos{\theta_{x}}\end{pmatrix}$
		$\displaystyle=\begin{pmatrix}\cos{\theta_{z}}\cos{\theta_{y}}&\cos{\theta_{z}}\sin{\theta_{y}}\sin{\theta_{x}}-\sin{\theta_{z}}\cos{\theta_{x}}&\cos{\theta_{z}}\sin{\theta_{y}}\cos{\theta_{x}}-\sin{\theta_{z}}\sin{\theta_{x}}\\ \sin{\theta_{z}}\cos{\theta_{y}}&\sin{\theta_{z}}\sin{\theta_{y}}\sin{\theta_{x}}+\cos{\theta_{z}}\cos{\theta_{x}}&\sin{\theta_{z}}\sin{\theta_{y}}\cos{\theta_{x}}-\cos{\theta_{z}}\sin{\theta_{x}}\\ -\sin{\theta_{y}}&\cos{\theta_{y}}\sin{\theta_{x}}&\cos{\theta_{y}}\cos{\theta_{x}}\end{pmatrix}$
		$\displaystyle=\begin{pmatrix}r_{11}&r_{12}&r_{13}\\ r_{21}&r_{22}&r_{23}\\ r_{31}&r_{32}&r_{33}\end{pmatrix}$

Then we obtain the Euler angles $\theta_{x},\theta_{y}\text{ and }\theta_{z}$ by

	$\displaystyle\theta_{x}$	$\displaystyle=atan2(r_{32},r_{33})$		(6)
	$\displaystyle\theta_{y}$	$\displaystyle=atan2(-r_{31},\sqrt{r_{32}^{2}+r_{33}^{2}})$
	$\displaystyle\theta_{z}$	$\displaystyle=atan2(r_{21},r_{11})$

In conclusion, we have the following two equivalent ways to define 3D quasiconformal representation for diffeomorphism:

1. 1.

   The upper triangular part of the matrix $P(p)$. $P(p)$ is a symmetric positive definite matrix, therefore the upper triangular part can encode the information of the matrix;

2. 2.

   The eigenvalues and the Euler angles $(a,b,c,\theta_{x},\theta_{y},\theta_{z})$. The Euler angles encode the information of the eigenvectors, therefore they include the equivalent information of $P(p)$.

Note that both representations include 6 numbers for each point in the domain. The first representation facilitates various computational operations, as the angles are not involved during the computation, while the second representation provides a much clearer geometric meaning. Since the two representations are equivalent and can be easily converted to each other, we will use only the first representation, i.e., the upper triangular part of $P(p)$, in the following sections without further mention.

By Equations 4, we have

$$\textbf{J}_{f}(p)^{T}\textbf{J}_{f}(p)=W\Sigma\Sigma W^{-1}=W\begin{pmatrix}a^{2}&0&0\\ 0&b^{2}&0\\ 0&0&c^{2}\end{pmatrix}W^{-1}$$

(7)

$$\det(\textbf{J}_{f}(p))=\det(P(p))=\det(\Sigma)=abc$$

(8)

Recall the assumption that $f$ is a diffeomorphism, implying $\textbf{J}_{f}(p)$ is invertible. Right multiplying by $\textbf{J}_{f}(p)^{-1}$ on both side of the above equation, we have

$$\textbf{J}_{f}(p)^{T}=W\begin{pmatrix}a^{2}&0&0\\ 0&b^{2}&0\\ 0&0&c^{2}\end{pmatrix}W^{-1}\frac{1}{\det(\textbf{J}_{f}(p))}C=W\begin{pmatrix}\frac{a}{bc}&0&0\\ 0&\frac{b}{ac}&0\\ 0&0&\frac{c}{ab}\end{pmatrix}W^{-1}C$$

(9)

where $C$ is the adjugate matrix of $\textbf{J}_{f}(p)$. According to the observation mentioned in Equation 2, $C$ has the form

$$C=\begin{pmatrix}|&|&|\\ \nabla v\times\nabla w&\nabla w\times\nabla u&\nabla u\times\nabla v\\ |&|&|\end{pmatrix}$$

(10)

where $\nabla u$, $\nabla v$ and $\nabla w$ are the columns of $\textbf{J}_{f}(p)^{T}$.

Moving the terms before the matrix C on the right-hand side of Equation 9 to the left, we have

$$W\begin{pmatrix}\frac{bc}{a}&0&0\\ 0&\frac{ac}{b}&0\\ 0&0&\frac{ab}{c}\end{pmatrix}W^{-1}\textbf{J}_{f}(p)^{T}=C$$

(11)

Let $\mathcal{A}$ be a matrix defined as

$$\mathcal{A}=W\begin{pmatrix}\frac{bc}{a}&0&0\\ 0&\frac{ac}{b}&0\\ 0&0&\frac{ab}{c}\end{pmatrix}W^{-1}$$

(12)

We can rewrite Equation 11 in the following form

$$\mathcal{A}\begin{pmatrix}|&|&|\\ \nabla u&\nabla v&\nabla w\\ |&|&|\end{pmatrix}=\begin{pmatrix}|&|&|\\ \nabla v\times\nabla w&\nabla w\times\nabla u&\nabla u\times\nabla v\\ |&|&|\end{pmatrix}$$

(13)

For $\nabla u$, we have

$$\nabla\cdot\mathcal{A}\nabla u=\nabla\cdot(\nabla v\times\nabla w)=\nabla w\cdot(\nabla\times\nabla v)-\nabla v\cdot(\nabla\times\nabla w)=0$$

(14)

since $\nabla\times\nabla$ results in zero vector. This also holds for $\nabla v$ and $\nabla w$.

Therefore, we compute divergence on both sides of Equation 13 to get the following equations

	$\displaystyle\nabla\cdot\mathcal{A}\nabla u$	$\displaystyle=0$		(15)
	$\displaystyle\nabla\cdot\mathcal{A}\nabla v$	$\displaystyle=0$
	$\displaystyle\nabla\cdot\mathcal{A}\nabla w$	$\displaystyle=0$

Furthermore, we define $\nabla\cdot\mathcal{A}\nabla$ as the 3D generalized Laplacian operator.

In the 2D case, we typically triangularize the domain to obtain finite elements and then define a piecewise linear function on that domain to approximate a smooth mapping. In the 3D case, we need to define the piecewise linear function on tetrahedrons serving as the elements.

Consider a 3D simply connected domain $\Omega_{1}$ discretized into a tetrahedral mesh containing vertices $\mathcal{V}=\{p_{1},p_{2},\ldots,p_{n}\}$ and tetrahedrons represented by the set of index sets $\mathcal{F}=\{(i_{1},i_{2},i_{3},i_{4})\mid i_{1},i_{2},i_{3},i_{4}\in\mathbb{N}\cap[1,n],\text{ and }i_{1},i_{2},i_{3},i_{4}\text{ are distinct}\}$. Also, the order of the vertices in a tetrahedron $T=(i_{1},i_{2},i_{3},i_{4})$ should satisfy a condition: the determinant of the matrix consisting of the three vectors $(p_{i_{2}}-p_{i_{1}})$, $(p_{i_{3}}-p_{i_{1}})$, and $(p_{i_{4}}-p_{i_{1}})$ must be positive, or equivalently, $[(p_{i_{2}}-p_{i_{1}})\times(p_{i_{3}}-p_{i_{1}})]\cdot(p_{i_{4}}-p_{i_{1}})>0$.

In this setting, the mapping $f:\Omega_{1}\rightarrow\Omega_{2}$ is piecewise linear.

To compute the 3D quasiconformal representation $\vec{q}$ for $f$, we need to compute $\vec{q}$ at each point in $\Omega$. Since $f$ is piecewise linear, the Jacobian matrices $\textbf{J}_{f}$ at each point in the same tetrahedron are equal to each other, resulting in the same polar decomposition and 3D quasiconformal representation. Therefore, we simply need to compute the Jacobian matrices at each tetrahedron $T$, denoted as $\textbf{J}_{f}(T)$, and use it to get $\vec{q}$. The linear function defined on each tetrahedron $T$ can be written as

$$f|_{T}(x,y,z)=\textbf{J}_{f}(T)\cdot\begin{pmatrix}x\\ y\\ z\end{pmatrix}+\begin{pmatrix}\vartheta_{T}\\ \iota_{T}\\ \varrho_{T}\end{pmatrix}$$

(16)

We can then compute the polar decomposition of $\textbf{J}_{f}(T)$ using Equations 4, which subsequently provides the 3D quasiconformal representation $\vec{q}$ for $T$.

To reconstruct the mapping $f$, we need the matrix $\mathcal{A}$, which contains the 3D representation information. Once given the representation $\vec{q}$ at a point $p$, $\mathcal{A}$ at $p$ can be computed according to Equation 12. As every point in $T$ shares the same dilation, $\mathcal{A}$ is a piecewise-constant matrix-valued function, where $\mathcal{A}$ takes a constant matrix value within each $T$, which we denote as $\mathcal{A}_{T}$.

To construct the discrete version of Equations 15, we should first derive the discrete gradient operator. For a certain tetrahedron $T=(i_{1},i_{2},i_{3},i_{4})$, let $p_{I}=(x_{I},y_{I},z_{I}),s_{I}=f(p_{I})=(u_{I},v_{I},w_{I})$, where $I\in\{i_{1},i_{2},i_{3},i_{4}\}$ . Then we have the equation

$$\textbf{J}_{f}(T)X=Y$$

(17)

where

$$X=\begin{pmatrix}|&|&|\\ p_{i_{2}}-p_{i_{1}}&p_{i_{3}}-p_{i_{1}}&p_{i_{4}}-p_{i_{1}}\\ |&|&|\end{pmatrix}\text{ and }Y=\begin{pmatrix}|&|&|\\ s_{i_{2}}-s_{i_{1}}&s_{i_{3}}-s_{i_{1}}&s_{i_{4}}-s_{i_{1}}\\ |&|&|\end{pmatrix}$$

(18)

Then moving the matrix $X$ to the right, we have

$$\textbf{J}_{f}(T)=YX^{-1}$$

(19)

For simplicity, we index the values in $X^{-1}$ by the form $\chi_{jk}$, where $j\text{ and }k$ are the row and column of the indexed entry.

We define

$\displaystyle\begin{gathered}A_{T}^{i_{1}}=-(\chi_{11}+\chi_{21}+\chi_{31}),\ A_{T}^{i_{2}}=\chi_{11},\ A_{T}^{i_{3}}=\chi_{21}\text{ and }A_{T}^{i_{4}}=\chi_{31}\\ B_{T}^{i_{1}}=-(\chi_{12}+\chi_{22}+\chi_{32}),\ B_{T}^{i_{2}}=\chi_{12},\ B_{T}^{i_{3}}=\chi_{22}\text{ and }B_{T}^{i_{4}}=\chi_{32}\\ C_{T}^{i_{1}}=-(\chi_{13}+\chi_{23}+\chi_{33}),\ C_{T}^{i_{2}}=\chi_{13},\ C_{T}^{i_{3}}=\chi_{23}\text{ and }C_{T}^{i_{4}}=\chi_{33}\\ \partial_{x}=(A_{T}^{i_{1}},A_{T}^{i_{2}},A_{T}^{i_{3}},A_{T}^{i_{4}})^{T},\partial_{y}=(B_{T}^{i_{1}},B_{T}^{i_{2}},B_{T}^{i_{3}},B_{T}^{i_{4}})^{T},\partial_{z}=(C_{T}^{i_{1}},C_{T}^{i_{2}},C_{T}^{i_{3}},C_{T}^{i_{4}})^{T}\\ u_{T}=(u_{i_{1}},u_{i_{2}},u_{i_{3}},u_{i_{4}})^{T},v_{T}=(v_{i_{1}},v_{i_{2}},v_{i_{3}},v_{i_{4}})^{T},w_{T}=(w_{i_{1}},w_{i_{2}},w_{i_{3}},w_{i_{4}})^{T}\end{gathered}$

(20)

Then we have the following result, from which we know that $\partial_{x}$, $\partial_{y}$, and $\partial_{z}$ compose the discrete gradient operator $\nabla_{T}$.

$$\textbf{J}_{f}(T)=\begin{pmatrix}\partial_{x}\cdot u_{T}&\partial_{y}\cdot u_{T}&\partial_{z}\cdot u_{T}\\ \partial_{x}\cdot v_{T}&\partial_{y}\cdot v_{T}&\partial_{z}\cdot v_{T}\\ \partial_{x}\cdot w_{T}&\partial_{y}\cdot w_{T}&\partial_{z}\cdot w_{T}\end{pmatrix}$$

(21)

The remaining part is the divergence operator. We first look at the definition of divergence. In a connected domain $\mathcal{W}$ , the divergence of a vector field F evaluated at a point $p\in\mathcal{W}$ is defined as

$$div\textbf{F}|_{p}=\lim_{|\mathcal{W}|\rightarrow 0}\frac{1}{|\mathcal{W}|}\oint_{\partial\mathcal{W}}\textbf{F}\cdot\hat{n}dS=\lim_{|\mathcal{W}|\rightarrow 0}\frac{1}{|\mathcal{W}|}\Phi(\textbf{F},\partial\mathcal{W})$$

(22)

where $|\mathcal{W}|$ denotes the volume of $\mathcal{W}$, $\partial\mathcal{W}$ is the boundary surface of $\mathcal{W}$, $\hat{n}$ is the outward unit normal to $\partial\mathcal{W}$, and $\Phi(\textbf{F},\partial\mathcal{W})$ is the flux of F across $\partial\mathcal{W}$.

In the discrete case, a vertex $p_{i}$ is surrounded by a set of tetrahedrons $\mathcal{N}(p_{i})$. Considering that our goal is not to develop a discrete divergence operator that can give precise divergence at every vertex, but to solve Equations 15 to get $u$, $v$, and $w$ such that the divergence of these functions is approximately equal to 0, we propose to replace the divergence operator by computing the flux in the discrete setting. The flux of F evaluated at $p$ can be written as follows

$$\Phi_{p_{i}}(\textbf{F},\partial\bigl{(}\cup\mathcal{N}(p_{i})\bigr{)})=\sum_{T\in\mathcal{N}(p_{i})}Area(T_{\setminus i})\cdot\textbf{F}(T)\cdot\hat{n}_{T_{\setminus i}}$$

(23)

here we use $T_{\setminus i}$ to denote the triangle opposite to $p_{i}$, $\hat{n}_{T_{\setminus i}}$ is the outward unit normal to the triangle $T_{\setminus i}$, and $\textbf{F}(T)$ represents the vector field defined on $T$. In our setting, $\mathbf{F}(T)$ is a piecewise-constant vector field, where $\mathbf{F}(T)$ takes a constant vector value within each $T$.

Considering that we shrink $T\in\mathcal{N}(p_{i})$ to $\mathcal{T}^{T}$ whose corresponding vertices $\{p_{i},p^{\mathcal{T}}_{2},p^{\mathcal{T}}_{3},p^{\mathcal{T}}_{4}\}$ satisfy $p^{\mathcal{T}}_{2}-p_{i}=\frac{1}{n}(p_{i_{2}}-p_{i_{1}})$, $p^{\mathcal{T}}_{3}-p_{i}=\frac{1}{n}(p_{i_{3}}-p_{i_{1}})$, and $p^{\mathcal{T}}_{4}-p_{i}=\frac{1}{n}(p_{i_{4}}-p_{i_{1}})$. Then the flux

	$\displaystyle\Phi_{p_{i}}(\textbf{F},\partial\bigl{(}\cup\{\mathcal{T}^{T}\mid T\in\mathcal{N}(p_{i})\}\bigr{)})$	$\displaystyle=\sum_{\mathcal{T}^{T},T\in\mathcal{N}(p_{i})}Area(\mathcal{T}^{T}_{\setminus i})\cdot\textbf{F}(\mathcal{T}^{T})\cdot\hat{n}_{\mathcal{T}^{T}_{\setminus i}}$
		$\displaystyle=\frac{1}{n^{2}}\sum_{T\in\mathcal{N}(p_{i})}Area(T_{\setminus i})\cdot\textbf{F}(T)\cdot\hat{n}_{T_{\setminus i}}$
		$\displaystyle=\frac{1}{n^{2}}\Phi_{p_{i}}(\textbf{F},\partial\bigl{(}\cup\mathcal{N}(p_{i})\bigr{)})$

since in the discrete setting, $\textbf{F}(\mathcal{T}^{T})=\textbf{F}(T)$, $\hat{n}_{\mathcal{T}^{T}_{\setminus i}}=\hat{n}_{T_{\setminus i}}$, and $Area(\mathcal{T}^{T}_{\setminus i})=\frac{1}{n^{2}}Area(T_{\setminus i})$. Therefore, if we require $\Phi_{p_{i}}(\textbf{F},\partial\bigl{(}\cup\mathcal{N}(p_{i})\bigr{)})=0$, then

	$\displaystyle\lim_{n\rightarrow\infty}\frac{1}{|\mathcal{W}(n)|}\Phi_{p_{i}}(\textbf{F},\partial\bigl{(}\cup\{\mathcal{T}^{T}\mid T\in\mathcal{N}(p_{i})\}\bigr{)})=$	$\displaystyle\lim_{n\rightarrow\infty}\frac{1}{|\mathcal{W}(n)|}\frac{1}{n^{2}}\Phi_{p_{i}}(\textbf{F},\partial\bigl{(}\cup\mathcal{N}(p_{i})\bigr{)})$
	$\displaystyle=$	$\displaystyle\lim_{n\rightarrow\infty}\frac{0}{n^{2}|\mathcal{W}(n)|}$
	$\displaystyle=$	$\displaystyle 0$

where $\mathcal{W}(n)=\cup\{\mathcal{T}^{T}\mid T\in\mathcal{N}(p_{i})\}$, and therefore, as $n\rightarrow\infty$, $|\mathcal{W}(n)|\rightarrow 0$. Although the above equation is not strictly equivalent to $div\textbf{F}|_{p}=0$, it motivates us to approximate $div\textbf{F}|_{p}=0$ using $\Phi_{p_{i}}(\textbf{F},\partial\bigl{(}\cup\mathcal{N}(p_{i})\bigr{)})=0$.

To solve the problem, we need to rewrite Equation 23 into a computable form. The quantity in the summation of Equation 23 can be evaluated using the property introduced in Equation 3. By comparing $X$ with the matrix $M$ in Equation 3, we notice that the three row vectors of $X^{-1}$ satisfy

	$\displaystyle(A_{T}^{i_{2}},B_{T}^{i_{2}},C_{T}^{i_{2}})^{T}$	$\displaystyle=\frac{1}{\det(X)}(p_{i_{3}}-p_{i_{1}})^{T}\times(p_{i_{4}}-p_{i_{1}})^{T}=-\frac{2Area(T_{\setminus i_{2}})}{6Vol(T)}\cdot\hat{n}_{T_{\setminus i_{2}}}$
	$\displaystyle(A_{T}^{i_{3}},B_{T}^{i_{3}},C_{T}^{i_{3}})^{T}$	$\displaystyle=\frac{1}{\det(X)}(p_{i_{4}}-p_{i_{1}})^{T}\times(p_{i_{2}}-p_{i_{1}})^{T}=-\frac{2Area(T_{\setminus i_{3}})}{6Vol(T)}\cdot\hat{n}_{T_{\setminus i_{3}}}$
	$\displaystyle(A_{T}^{i_{4}},B_{T}^{i_{4}},C_{T}^{i_{4}})^{T}$	$\displaystyle=\frac{1}{\det(X)}(p_{i_{2}}-p_{i_{1}})^{T}\times(p_{i_{3}}-p_{i_{1}})^{T}=-\frac{2Area(T_{\setminus i_{4}})}{6Vol(T)}\cdot\hat{n}_{T_{\setminus i_{4}}}$

where $Vol(T)$ represents the volume of $T$ with respect to $\mathcal{V}$, and $\det(X)=6Vol(T)$.

For vertex $p_{i_{1}}$, the outward normal vector on $T_{\setminus i_{1}}$ is given by

	$\displaystyle(p_{i_{3}}-p_{i_{2}})^{T}\times(p_{i_{4}}-p_{i_{2}})^{T}=$	$\displaystyle[(p_{i_{3}}-p_{i_{1}})^{T}-(p_{i_{2}}-p_{i_{1}})^{T}]\times[(p_{i_{4}}-p_{i_{1}})^{T}-(p_{i_{2}}-p_{i_{1}})^{T}]$
	$\displaystyle=$	$\displaystyle(p_{i_{3}}-p_{i_{1}})^{T}\times(p_{i_{4}}-p_{i_{1}})^{T}-(p_{i_{3}}-p_{i_{1}})^{T}\times(p_{i_{2}}-p_{i_{1}})^{T}$
		$\displaystyle-(p_{i_{2}}-p_{i_{1}})^{T}\times(p_{i_{4}}-p_{i_{1}})^{T}+(p_{i_{2}}-p_{i_{1}})^{T}\times(p_{i_{2}}-p_{i_{1}})^{T}$
	$\displaystyle=$	$\displaystyle\det(X)[(A_{T}^{i_{2}},B_{T}^{i_{2}},C_{T}^{i_{2}})^{T}+(A_{T}^{i_{3}},B_{T}^{i_{3}},C_{T}^{i_{3}})^{T}+(A_{T}^{i_{4}},B_{T}^{i_{4}},C_{T}^{i_{4}})^{T}]$
	$\displaystyle=$	$\displaystyle\det(X)(A_{T}^{i_{2}}+A_{T}^{i_{3}}+A_{T}^{i_{4}},B_{T}^{i_{2}}+B_{T}^{i_{3}}+B_{T}^{i_{4}},C_{T}^{i_{2}}+C_{T}^{i_{3}}+C_{T}^{i_{4}})^{T}$
	$\displaystyle=$	$\displaystyle-\det(X)(A_{T}^{i_{1}},B_{T}^{i_{1}},C_{T}^{i_{1}})^{T}$

We also know that $(p_{i_{3}}-p_{i_{2}})^{T}\times(p_{i_{4}}-p_{i_{2}})^{T}=2Area(T_{\setminus i_{1}})\cdot\hat{n}_{T_{\setminus i_{1}}}$. The geometric meaning of $-(A_{T}^{i},B_{T}^{i},C_{T}^{i})^{T}$ is shown in Figure 2.

Together, for the term in the summation of flux evaluated at the four vertices of $T$, we have

$$Area(T_{\setminus i})\cdot\textbf{F}(T)\cdot\hat{n}_{T_{\setminus i}}=-3Vol(T)\cdot(A_{T}^{i},B_{T}^{i},C_{T}^{i})^{T}\cdot\textbf{F}(T),\quad i\in T$$