Version of August 28, 2024$,\,\,\,$

# Nonnegative moment coordinates on finite element geometries

Luca Dieci  School of Mathematics, Georgia Institute of Technology, Atlanta,
GA 30332 U.S.A<EMAIL_ADDRESS>,  Fabio V. Difonzo  Dipartimento di
Matematica, Università degli Studi di Bari Aldo Moro, Via E. Orabona 4, 70125
Bari, Italy<EMAIL_ADDRESS>and  N. Sukumar  Department of Civil and
Environmental Engineering, One Shields Avenue, University of California,
Davis, CA 95616, USA<EMAIL_ADDRESS>
###### Abstract.

In this paper, we introduce new generalized barycentric coordinates (coined as
moment coordinates) on nonconvex quadrilaterals and convex hexahedra with
planar faces. This work draws on recent advances in constructing interpolants
to describe the motion of the Filippov sliding vector field in nonsmooth
dynamical systems, in which nonnegative solutions of signed matrices based on
(partial) distances are studied. For a finite element with $n$ vertices
(nodes) in $\mathbb{R}^{2}$, the constant and linear reproducing conditions
are supplemented with additional linear moment equations to set up a linear
system of equations of full rank $n$, whose solution results in the
nonnegative shape functions. On a simple (convex or nonconvex) quadrilateral,
moment coordinates using signed distances are identical to mean value
coordinates. For signed weights that are based on the product of distances to
edges that are incident to a vertex and their edge lengths, we recover
Wachspress coordinates on a convex quadrilateral. Moment coordinates are also
constructed on a convex hexahedra with planar faces. We present proofs in
support of the construction and plots of the shape functions that affirm its
properties.

###### Key words and phrases:

generalized barycentric coordinates; mean value coordinates; Wachspress
coordinates; quadrilateral; hexahedron; Filippov vector field

###### Key words and phrases:

generalized barycentric coordinates; mean value coordinates; Wachspress
coordinates; quadrilateral; hexahedron; Filippov vector field

###### 1991 Mathematics Subject Classification:

52A07, 52B55

##  1\. Introduction

Generalized barycentric coordinates (GBCs) [1, 2] are used to interpolate data
that are prescribed at the vertices of a polytope. These coordinates (shape
functions) extend the well-known barycentric coordinates on simplices to
arbitrary polytopes in $\mathbb{R}^{d}$ with more than $d+1$ vertices. GBCs
have found applications in computer graphics, computational mechanics, and the
applied sciences [3]. Let $P\subset\mathbb{R}^{d}$ be a polytope (convex or
nonconvex) with $n$ vertices $\\{v_{i}\\}_{i=1}^{n}$. A point in $P$ is
denoted by $p\in\mathbb{R}^{d}$. In its most general definition, generalized
barycentric coordinates, $\varphi(p):=\\{\varphi_{i}\\}_{i=1}^{n}$, must
satisfy the following relations:

(1.1)
$\sum_{i=1}^{n}\varphi_{i}(p)=1,\quad\sum_{i=1}^{n}\varphi_{i}(p)v_{i}=p,\quad\varphi_{i}(v_{j})=\delta_{ij},\quad\varphi_{i}\geq
0,$

where $\delta_{ij}$ is the Kronecker-delta and (1.1) ensures that the
generalized barycentric interpolant has linear precision. In addition, it is
also desirable that $\varphi$ are smooth in the interior of $P$.

Wachspress [4] was the first to introduce nonnegative rational basis that were
valid for convex polygons; Warren [5] extended it for convex polyhedra.
Floater [6] proposed the mean value coordinates, which have a simple
expression that is amenable to numerical computations, and is the most popular
and widely used set of GBCs. These coordinates are valid for simple (convex
and nonconvex) and nested polygons, but the coordinates can become negative on
nonconvex polygons [7]. For applications such as shape deformation in computer
graphics, data approximation over meshes that preserve the maximum principle
and mitigate the Runge phenomenon, and to ensure positivity of the mass matrix
entries in an explicit Galerkin method for elastodynamic simulations, it is
beneficial that the nonnegativity condition on $\varphi$ is met; however, over
the years, many coordinates have been proposed that can become negative in $P$
[2]. On convex polytopes, coordinates that are nonnegative are endowed with
the facet-reducing property on its faces [8]. The only coordinates that are
nonnegative and can be $C^{1}$ smooth in $P$ are harmonic coordinates [9],
maximum-entropy coordinates [10, 11], local barycentric coordinates [12, 13],
blended barycentric coordinates [14] and iterative coordinates [15]. Harmonic
coordinates require the numerical solution of the Laplace equation with
piecewise affine boundary conditions, whereas those that follow need numerical
computations (not known analytically) and/or use a partitioning of the
polytope.

The construction of generalized barycentric coordinates that satisfy (1.1) and
are given by a simple analytical expression remains an open problem. In this
contribution, we propose an initial attempt to fill this void for common
finite element geometries. To this end, we supplement the $d+1$ constraints in
(1.1) by additional moment constraints so that the resulting linear system has
a unique nonnegative solution, which is the desired $\varphi$. We refer to the
GBCs so constructed as either moment coordinates (MC) or moment shape
functions. Herein, we develop these coordinates for finite element shapes:
quadrilaterals in $\mathbb{R}^{2}$ and hexahedra in $\mathbb{R}^{3}$. These
geometries can be simple quadrilaterals (convex or nonconvex) and convex
hexahedra. This idea of adding moment equality conditions to (1.1) traces back
to sliding mode control in piecewise smooth dynamical systems (see, for
completeness, [16] for the quadrilateral in $\mathbb{R}^{2}$ and [17] for its
extension to the hexahedron in $\mathbb{R}^{3}$), when one looks for Filippov
solutions in presence of dynamics along discontinuity manifolds of codimension
2 or higher. The connections of the moment approach to mean value coordinates
are discussed in [16].

The structure of the remainder of this paper follows. In Section 2, the
formulation to derive moment coordinates is presented: quadrilaterals are
treated in Sections 2.1 and 2.2, and piecewise affine interpolation in the
one-dimensional setting is considered in Section 2.3. Extension of moment
coordinates to hexahedra is presented is Section 2.4. Numerical results in
Section 3 include closed-form expressions and plots of the moment coordinates
on finite element geometries. We summarize our main finding with some
perspectives on future work in Section 4.

##  2\. Formulation

We first describe the approach to obtain analytical (using symbolic
computations) expressions of the moment coordinates on simple quadrilaterals
in Section 2.1 (recover mean value coordinates) and Section 2.2 (recover
Wachspress coordinates). Then we consider one-dimensional piecewise affine
interpolation in Section 2.3 and moment coordinates on convex hexahedra in
Section 2.4. Our main results provide invertibility for matrices with
prescribed sign pattern; however, such matrices elude classical tools for
proving invertibility (for example, see [18, 19]) and would require specific
analysis. We also point out that the sign pattern is identical to the signs of
the hourglass nodal vectors (spurious zero-energy modes of the stiffness
matrix for the Poisson equation) for the four-node quadrilateral and the
eight-node hexahedral element [20].

###  2.1. Four-node (convex or nonconvex) quadrilateral

$v_{1}$$v_{2}$$v_{3}$$v_{4}$$p$$d_{1}$$d_{2}$$d_{3}$$d_{4}$$h_{1}$$h_{2}$$h_{3}$$h_{4}$$n_{1}$$n_{2}$$n_{3}$$n_{4}$
(a)

$v_{1}$$v_{2}$$v_{3}$$v_{4}$$p$$d_{1}$$d_{2}$$d_{3}$$d_{4}$ (b)

Figure 1. Simple quadrilaterals. (a) Convex and (b) Nonconvex.
$\\{v_{i}\\}_{i=1}^{4}$ are the vertices of the quadrilateral $Q$ and $p$ is a
generic point in $Q$. Vertex ordering is assumed to be cyclic, so that
$v_{5}=v_{1}$ and $v_{0}=v_{4}$.

Let $\\{v_{i}\\}_{i=1}^{4}\subseteq\mathbb{R}^{2}$ be affinely independent
vertices of a simple quadrilateral (see Figure 1). The set of all possible
barycentric coordinates for any $p\in Q$,
$Q\mathrel{\mathop{:}}=\textrm{conv}\,\\{v_{i}\\}$, is the set of all the
nonnegative solutions $\varphi(p)$ to the full-rank, underdetermined system

(2.1) $\begin{bmatrix}\mathds{1}^{\top}\\\
V\end{bmatrix}\varphi(p)=\begin{bmatrix}1\\\ p\end{bmatrix}.$

As a consequence of the affine independence, such a system has a one-
dimensional kernel, so that there exists a nonzero $g\in\mathbb{R}^{4}$ and
$\langle g\rangle=\ker\begin{bmatrix}\mathds{1}^{\top}\\\ V\end{bmatrix}$. It
is straightforward to see that $g$ has a specific sign pattern, which can be
computed by

$g=\begin{bmatrix}\tau(v_{4})\\\ -1\end{bmatrix},$

where $\tau(v_{4})$ represents the triangular barycentric coordinate vector of
$v_{4}$ with respect to the triangle $v_{1},v_{2},v_{3}$. Thus, it turns out
that

$\mathrm{sgn}(g)=\begin{bmatrix}+\\\ -\\\ +\\\ -\end{bmatrix}.$

Therefore, if we want to select a specific set of barycentric coordinates for
each $p\in Q$, one approach is to suitably append an additional condition to
(2.1) of the form

$d(p)^{\top}\varphi(p)=0,\quad d(p)\in\mathbb{R}^{4},$

so as to make it nonsingular and preserve convexity for its unique solution.
Of course, if the vector $d(p)$ is not orthogonal to $g$ then the resulting
system would be nonsingular. The easiest way to accomplish this is to select
$d(p)$ such that $\mathrm{sgn}(d)=\mathrm{sgn}(g)$. Now, for convexity to be
preserved, we would at least guarantee it holds on $\partial Q$. To this end,
let $p\in v_{1}v_{2}$, so that we expect

$\varphi(p)=\begin{bmatrix}1-\alpha\\\ \alpha\\\ 0\\\
0\end{bmatrix},\quad\alpha\in[0,1],$

to be the unique solution to

(2.2) $\begin{bmatrix}\mathds{1}^{\top}\\\ V\\\
d(p)^{\top}\end{bmatrix}\varphi(p)=\begin{bmatrix}1\\\ p\\\ 0\end{bmatrix}.$

But since $\alpha=\frac{\|p-v_{1}\|}{\|v_{2}-v_{1}\|}$, we get

$\frac{\|p-v_{2}\|}{\|v_{2}-v_{1}\|}d_{1}+\frac{\|p-v_{1}\|}{\|v_{2}-v_{1}\|}d_{2}=0,$

from which, if
$d_{1}\mathrel{\mathop{:}}=\|p-v_{1}\|,d_{2}\mathrel{\mathop{:}}=-\|p-v_{2}\|$,
the above holds true (see Figures 1(a) and 1(b)). Repeating the same argument
for the remaining edges, we see that defining

(2.3)
$d(p)\mathrel{\mathop{:}}=\begin{bmatrix}d_{1}(p)&-d_{2}(p)&d_{3}(p)&-d_{4}(p)\end{bmatrix},\quad
d_{i}(p)\mathrel{\mathop{:}}=\|p-v_{i}\|,\,\,i=1,2,3,4,$

provides (2.2) to be nonsingular and its unique solution to be convex on
$\partial Q$. Proving that it is also feasible in the interior of $Q$ can be
found in [16]. An alternative proof for the convex quadrilateral case follows.

######  Proposition 2.1.

For any $p\in Q$ there exist $\alpha,\beta,\gamma\in[0,1]$ such that

$p=(1-\gamma)a+\gamma b,\quad a\mathrel{\mathop{:}}=(1-\alpha)v_{4}+\alpha
v_{1},\,\,b\mathrel{\mathop{:}}=(1-\beta)v_{2}+\beta v_{3}$

and

(2.4)
$d(p)^{\top}\varphi(p)=0,\quad\varphi(p)\mathrel{\mathop{:}}=\begin{bmatrix}(1-\gamma)\alpha&\gamma(1-\beta)&\gamma\beta&(1-\gamma)(1-\alpha)\end{bmatrix}^{\top},$

where $d(p)$ is as in (2.3).

###### Proof.

Let us assume that $p=0$, and without loss of generality, $d_{i}=1$ for all
$i=1,2,3,4$. Therefore, (2.4) reads

$\alpha-\gamma\alpha+\gamma\beta=\frac{1}{2}.$

The function
$f(\alpha,\beta,\gamma)\mathrel{\mathop{:}}=\alpha-\gamma\alpha+\gamma\beta$
is such that $\nabla
f(\alpha,\beta,\gamma)=\begin{bmatrix}1-\gamma&\gamma&\beta-\alpha\end{bmatrix}$,
so that $\nabla f\neq 0$ on $\mathbb{R}^{3}$. Moreover, on the boundary of
$[0,1]^{3}$ it attains the following values:

$\displaystyle f(0,\beta,\gamma)$ $\displaystyle=\gamma\beta,$ $\displaystyle
f(1,\beta,\gamma)$ $\displaystyle=1-\gamma(1-\beta),$ $\displaystyle
f(\alpha,0,\gamma)$ $\displaystyle=\alpha(1-\gamma),$ $\displaystyle
f(\alpha,1,\gamma)$ $\displaystyle=\alpha+\gamma(1-\alpha),$ $\displaystyle
f(\alpha,\beta,0)$ $\displaystyle=\alpha,$ $\displaystyle f(\alpha,\beta,1)$
$\displaystyle=\beta,$

and so has minimum $0$ and maximum $1$. By the Intermediate Value Theorem
there exist $\alpha,\beta,\gamma\in[0,1]$ such that
$f(\alpha,\beta,\gamma)=\frac{1}{2}$, which proves the claim. ∎

Even though moment coordinates are identical to mean value coordinates on a
quadrilateral, solving the linear system (2.2) is computationally attractive
to obtain $\varphi(p)$ for any point $p\in Q$. This is so, since the local
form of computing the mean value coordinates cannot be used if $p\in\partial
Q$ [1].

An alternative expression for moment barycentric coordinates on quadrilaterals
can be given in terms of triangular barycentric coordinates, which can be seen
as a trivial extension of barycentric coordinates relative to one of the
triangles when a quadrilateral is split by its diagonals, with a zero
component in the position of the removed vertex. For example, leaving $v_{4}$
out and considering the triangle $v_{1},v_{2},v_{3}$, then the corresponding
triangular barycentric coordinates $\tau^{1}(p)$ are computed for each $p\in
Q$ as

$\begin{bmatrix}1&1&1&1\\\ v_{1}&v_{2}&v_{3}&v_{4}\\\
1&0&0&0\end{bmatrix}\tau^{1}(p)=\begin{bmatrix}1\\\ p\\\ 0\end{bmatrix}.$

Applying Cramer’s rule to (2.2) for computing $\varphi_{1}(p)$ provides

$\varphi_{1}(p)=\frac{\det\begin{bmatrix}1&1&1&1\\\ p&v_{2}&v_{3}&v_{4}\\\
0&-d_{2}&d_{3}&-d_{4}\end{bmatrix}}{d(p)^{\top}\nu},$

where $\nu=\begin{bmatrix}A_{1}\\\ -A_{2}\\\ A_{3}\\\ -A_{4}\end{bmatrix}$,
$A_{i}\mathrel{\mathop{:}}=\lvert\mathrm{area}(\textrm{conv}\,\\{v_{j}\,:\,j\neq
i\\})\rvert$, and $\langle\nu\rangle=\ker\begin{bmatrix}\mathds{1}^{\top}\\\
V\end{bmatrix}$. Thus, observing that
$\nu^{1}(p)\mathrel{\mathop{:}}=A_{1}\begin{bmatrix}-1\\\
\tau^{1}(p)\end{bmatrix}$ is such that

$\langle\nu^{1}(p)\rangle=\ker\begin{bmatrix}1&1&1&1\\\
p&v_{2}&v_{3}&v_{4}\end{bmatrix},$

we have

$\displaystyle\varphi_{1}(p)$
$\displaystyle=\frac{\begin{bmatrix}0&-d_{2}&d_{3}&-d_{4}\end{bmatrix}\nu^{1}(p)}{d(p)^{\top}\nu}$
$\displaystyle=A_{1}\frac{\begin{bmatrix}0&-d_{2}&d_{3}&-d_{4}\end{bmatrix}\begin{bmatrix}-1\\\
\tau^{1}(p)\end{bmatrix}}{d(p)^{\top}\nu}$
$\displaystyle=A_{1}\frac{\begin{bmatrix}d_{1}&-d_{2}&d_{3}&-d_{4}\end{bmatrix}\begin{bmatrix}0\\\
\tau^{1}(p)\end{bmatrix}}{d(p)^{\top}\nu}$
$\displaystyle=A_{1}\frac{d(p)^{\top}\tau^{1}(p)}{d(p)^{\top}\nu}.$

Analogous arguments hold for the remaining components. Thus, in general, the
following holds:

$\varphi_{i}(p)=(-1)^{d+1+i}A_{i}\frac{d(p)^{\top}\tau^{i}(p)}{d(p)^{\top}\nu}.$

###  2.2. Convex quadrilaterals

In this section we prove that on a convex quadrilateral $Q$ (see Figure 1(a)),
akin to mean values coordinates, it is possible to compute Wachspress
coordinates by solving a linear system. Let us recall from [21] that
Wachspress coordinates can be expressed as

$\varphi(p)=\frac{w_{i}(p)}{\sum_{j=1}^{4}w_{j}(p)},$

where

(2.5)
$w_{i}(p)=\frac{A(v_{i-1},v_{i},v_{i+1})}{A(p,v_{i-1},v_{i})A(p,v_{i},v_{i+1})},$

and $A(a,b,c)$ stands for the signed area of the triangle with vertices
$a,b,c$ given by

$A(a,b,c)\mathrel{\mathop{:}}=\frac{1}{2}\det\begin{bmatrix}a_{1}&b_{1}&c_{1}\\\
a_{2}&b_{2}&c_{2}\\\ 1&1&1\end{bmatrix}.$

As is well known (see [1]), Wachspress coordinates can also, and more
conveniently, be expressed in terms of distances from the edges. More
precisely, for $i=1,2,3,4$ and with cyclic convention, letting $n_{i}$ be the
outer normal to the edge $v_{i}v_{i+1}$ and $h_{i}(p)$ be the distance from
$p$ to the edge $v_{i}v_{i+1}$, we define the weights

$\tilde{w}_{i}(p)\mathrel{\mathop{:}}=\frac{n_{i-1}\times
n_{i}}{h_{i-1}(p)h_{i}(p)},$

where, for general $x,y\in\mathbb{R}^{2}$, we define

$x\times y\mathrel{\mathop{:}}=\begin{vmatrix}x_{1}&x_{2}\\\
y_{1}&y_{2}\end{vmatrix}.$

It then follows that
$A(v_{i-1},v_{i},v_{i+1})=\frac{1}{2}l_{i,i+1}l_{i,i-1}n_{i-1}\times n_{i}$,
where we let $l_{i,i+1}\mathrel{\mathop{:}}=\|v_{i}-v_{i+1}\|$, for
$i=1,2,3,4$. It can be seen that $w_{i}(p)=2\tilde{w}_{i}(p)$. Now, letting

$\rho_{i}(p)\mathrel{\mathop{:}}=l_{i,i-1}l_{i,i+1}h_{i-1}(p)h_{i}(p),$

and using planar geometry, it is readily shown that

$\rho_{1}(p)\varphi_{1}(p)-\rho_{2}(p)\varphi_{2}(p)+\rho_{3}(p)\varphi_{3}(p)-\rho_{4}(p)\varphi_{4}(p)=0.$

Thus, letting

$\rho(p)\mathrel{\mathop{:}}=\begin{bmatrix}\rho_{1}(p)&-\rho_{2}(p)&\rho_{3}(p)&-\rho_{4}(p)\end{bmatrix}^{\top},$

it follows that $\rho(p)^{\top}v\neq 0$ because of the sign pattern, and so
Wachspress coordinates are the unique solution to the linear system

(2.6) $\begin{bmatrix}V\\\ \mathds{1}^{\top}\\\
\rho(p)^{\top}\end{bmatrix}\varphi(p)=\begin{bmatrix}p\\\ 1\\\
0\end{bmatrix}.$

We point out that $p\in Q$, $p\in v_{i}v_{i+1}$ if and only if $h_{i}(p)=0$.
Therefore, $p\in v_{i-1}v_{i}\cup v_{i}v_{i+1}$ if and only if
$\rho_{i}(p)=0$.

###  2.3. One-dimensional case

Moment condition can be leveraged to recover piecewise linear basis functions
on an interval in one dimension. Let us start with the case of three points in
$[0,1]$, say $x_{1}=0,x_{2}\in(0,1),x_{3}=1$. Letting $x\in[0,1]$, we then
have to consider the underdetermined linear system

(2.7) $\begin{bmatrix}1&1&1\\\
x_{1}&x_{2}&x_{3}\end{bmatrix}\varphi=\begin{bmatrix}1\\\ x\end{bmatrix}.$

$p$$x_{1}$$x_{2}$$x_{3}$$x$$A$$B$ Figure 2. Three nodes in one dimension.

The idea is to regularize the above system so that the point $(x,0)$ belongs
to the edge $AB$ (see Figure 2), where $A(x_{0},x-x_{1}),B(x_{1},x-x_{2}))$,
while ensuring that (2.7) is invertible so that its unique solution is
componentwise nonnegative. One way to accomplish this is to append to (2.7) an
additional row using moment regularization. Specifically, we propose to
regularize (2.7) as

$\begin{bmatrix}1&1&1\\\ x_{1}-x&x_{2}-x&x_{3}-x\\\
|x_{1}-x|&-|x_{2}-x|&|x_{3}-x|\end{bmatrix}\varphi=\begin{bmatrix}1\\\ 0\\\
0\end{bmatrix}.$

The system above is easily seen to be invertible, and since
$\begin{bmatrix}0\\\ 0\end{bmatrix}$ is in the convex hull of the points
$\begin{bmatrix}x_{1}-x\\\ |x_{1}-x|\end{bmatrix},\begin{bmatrix}x_{2}-x\\\
-|x_{2}-x|\end{bmatrix},\begin{bmatrix}x_{3}-x\\\ |x_{3}-x|\end{bmatrix}$,
then its unique solution $\varphi(x)$ is necessarily feasible.

In the following we extend such an argument to any number $n$ of points
$\\{v_{i}\\}_{i=1}^{n}$ in $[0,1]$. We recall a handy result from [17].

######  Lemma 2.2.

Let $A\in\mathbb{R}^{n\times m}$, $n<m$, be full rank, and let
$b\in\mathbb{R}^{n}$. Consider the system

(2.8) $Ax=b\ ,$

and let $d\in\mathbb{R}^{m}$ be a nonzero vector.

If there exist $x$ and $y$ solutions to (2.8), such that

$d^{\top}x=\xi\ ,\quad\text{and}\quad d^{\top}y=\eta\ ,$

with $\xi\neq\eta$, then $\begin{bmatrix}A\\\ d^{\top}\end{bmatrix}$ has rank
$n+1$.

Let $n\geq 3$,
$V\mathrel{\mathop{:}}=\begin{bmatrix}x_{1}&\cdots&x_{n}\end{bmatrix}$, with
$x_{i}\in[0,1]$ for $i=1,\ldots,n$, not necessarily evenly spaced. Then, given
$p\in[0,1]$, we seek for a feasible solution $\varphi(p)$ of the
underdetermined system

(2.9) $\begin{bmatrix}1&\cdots&1&\cdots&1\\\
x_{1}&\cdots&x_{i}&\cdots&x_{n}\\\
|x_{1}-x|&\cdots&(-1)^{i+1}|x_{i}-x|&\cdots&(-1)^{n+1}|x_{n}-x|\end{bmatrix}\varphi=\begin{bmatrix}1\\\
x\\\ 0\end{bmatrix},$

which has full rank $3$.

Without loss of generality, let us assume $p\in[x_{1},x_{2}]$; other cases can
be handled analogously. Let us note that the last row above, which may be
referred to as _moment regularization row_ , is always compatible with the
expected solution to (2.9), which has the form

$\varphi_{i}(p)=\begin{cases}(1-\alpha),&i=1,\\\ \alpha,&i=2,\\\ 0,&i\neq
1,2.\end{cases}$

The case $n=4$ is the classical moment system, which is known to be
invertible, providing a unique feasible solution. For $n\geq 5$, let us define

$\displaystyle d_{i}(p)$
$\displaystyle\mathrel{\mathop{:}}=(-1)^{i+1}|x_{i}-p|,\quad i=1,\ldots,n,$
$\displaystyle\Delta$
$\displaystyle\mathrel{\mathop{:}}=\begin{bmatrix}0&0&1&1&\cdots&0\\\
\vdots&\vdots&&\ddots&\ddots&\vdots\\\ 0&0&\cdots&\cdots&1&1\end{bmatrix}.$

We prove the following.

######  Proposition 2.3.

Let $p\in[0,1]$. For all $n\geq 4$ the system

$\begin{bmatrix}\mathds{1}^{\top}\\\ V-p\\\ d(p)^{\top}\\\
\Delta\end{bmatrix}\varphi(p)=\begin{bmatrix}1\\\ 0\\\ 0\\\
0_{n-3}\end{bmatrix}$

is invertible and its unique solution is feasible.

###### Proof.

We proceed by induction on $n$. For $n=4$

$\Delta=\begin{bmatrix}0&0&1&1\end{bmatrix}.$

Assume $p=x_{1}$. Let $\alpha\neq 0$ such that

$p=(1-\alpha)x_{2}+\alpha x_{3}.$

Therefore $0=\Delta\begin{bmatrix}1\\\ 0\\\ 0\\\
0\end{bmatrix}\neq\Delta\begin{bmatrix}0\\\ 1-\alpha\\\ \alpha\\\
0\end{bmatrix}=\alpha$, and Lemma 2.2 proves the claim. An analogous argument
hold if $p=x_{2}$. If $p\in(x_{1},x_{2})$, then let $\alpha,\beta\in(0,1)$
such that

$p=(1-\alpha)x_{1}+\alpha x_{2}=(1-\beta)x_{1}+\beta x_{4}.$

Again $0=\Delta\begin{bmatrix}1-\alpha\\\ \alpha\\\ 0\\\
0\end{bmatrix}\neq\Delta\begin{bmatrix}1-\beta\\\ 0\\\ 0\\\
\beta\end{bmatrix}=\beta$, and Lemma 2.2 proves the claim.
Let the result hold true for $n-1$, and let us define

$\delta^{\top}\mathrel{\mathop{:}}=\begin{bmatrix}0&0&\cdots&0&1&1\end{bmatrix}\in\mathbb{R}^{n}.$

Let $p\in(x_{1},x_{2})$. Then there exist $\alpha,\beta\in(0,1)$ such that

$p=(1-\alpha)x_{1}+\alpha x_{2}$

and

$p=\begin{cases}(1-\beta)x_{1}+\beta x_{n-1},&\textrm{if $n$ is odd},\\\
(1-\beta)x_{1}+\beta x_{n},&\textrm{if $n$ is even}.\end{cases}$

Therefore, assuming $n$ to be even,
$\varphi_{12}\mathrel{\mathop{:}}=\begin{bmatrix}1-\alpha\\\ \alpha\\\ 0\\\
\vdots\\\ 0\end{bmatrix},\
\varphi_{1n}\mathrel{\mathop{:}}=\begin{bmatrix}1-\beta\\\ 0\\\ \vdots\\\ 0\\\
\beta\end{bmatrix}$ are solutions to

$\begin{bmatrix}\mathds{1}^{\top}&1\\\ V_{1:n-1}-p&x_{n}-p\\\
d(p)^{\top}&-|x_{n}-x|\\\ \Delta&0\end{bmatrix}\varphi(p)=\begin{bmatrix}1\\\
0\\\ 0\\\ 0_{n-3}\end{bmatrix},$

and $\delta^{\top}\varphi_{12}=0\neq\beta=\delta^{\top}\varphi_{1n}$, so that
Lemma 2.2 proves the result. Similar arguments hold for the boundary cases
$p=x_{1}$ or $p=x_{2}$. We proceed analogously for $n$ odd. The claim is
proved. ∎

###  2.4. Eight-node hexahedron

We begin by first stating the main result from [17].

######  Theorem 2.4.

Let $v_{i}=\left[\begin{smallmatrix}v_{i}^{1}\\\ v_{i}^{2}\\\
v_{i}^{3}\end{smallmatrix}\right]$, $i=1,\dots,8$, be eight vectors in
$\mathbb{R}^{3}$, and consider the matrix $V\in\mathbb{R}^{3\times 8}$ given
by

(2.10)
$V\mathrel{\mathop{:}}=\begin{bmatrix}v_{1}&v_{2}&v_{3}&v_{4}&v_{5}&v_{6}&v_{7}&v_{8}\end{bmatrix}\
.$

Let $p\in H$, $H:=\textrm{conv}\,\\{v_{i}\\}$. Assume that the entries of
$V-p$ are nonzero and have the following signs:

(2.11) $\begin{bmatrix}+&+&+&+&-&-&-&-\\\ +&+&-&-&+&+&-&-\\\
+&-&-&+&+&-&-&+\end{bmatrix}.$

Let $\Delta(p)\in\mathbb{R}^{3\times 8}$ be the following matrix of signed
partial distances:

(2.12)
$\Delta(p)\mathrel{\mathop{:}}=\begin{bmatrix}\delta^{23}_{1}(p)&-\delta^{23}_{2}(p)&\delta^{23}_{3}(p)&-\delta^{23}_{4}(p)&\delta^{23}_{5}(p)&-\delta^{23}_{6}(p)&\delta^{23}_{7}(p)&-\delta^{23}_{8}(p)\\\
\delta^{13}_{1}(p)&-\delta^{13}_{2}(p)&-\delta^{13}_{3}(p)&\delta^{13}_{4}(p)&-\delta^{13}_{5}(p)&\delta^{13}_{6}(p)&\delta^{13}_{7}(p)&-\delta^{13}_{8}(p)\\\
\delta^{12}_{1}(p)&\delta^{12}_{2}(p)&-\delta^{12}_{3}(p)&-\delta^{12}_{4}(p)&-\delta^{12}_{5}(p)&-\delta^{12}_{6}(p)&\delta^{12}_{7}(p)&\delta^{12}_{8}(p)\end{bmatrix},$

where for each $i=1,\ldots,8$,

$\displaystyle\delta^{23}_{i}(p)$
$\displaystyle\mathrel{\mathop{:}}=\sqrt{(v_{i}^{2}-p_{2})^{2}+(v_{i}^{3}-p_{3})^{2}}\
,$ $\displaystyle\delta^{13}_{i}(p)$
$\displaystyle\mathrel{\mathop{:}}=\sqrt{(v_{i}^{1}-p_{1})^{2}+(v_{i}^{3}-p_{3})^{2}}\
,$ $\displaystyle\delta^{12}_{i}(p)$
$\displaystyle\mathrel{\mathop{:}}=\sqrt{(v_{i}^{1}-p_{1})^{2}+(v_{i}^{2}-p_{2})^{2}}\
.$

Finally, for each $i=1,\ldots,8$, let

$d(p)^{\top}\mathrel{\mathop{:}}=\begin{bmatrix}d_{1}(p)&-d_{2}(p)&d_{3}(p)&-d_{4}(p)&-d_{5}(p)&d_{6}(p)&-d_{7}(p)&d_{8}(p)\end{bmatrix},\quad
d_{i}(p)\mathrel{\mathop{:}}=\|v_{i}-p\|_{2},$

and let $\mathds{1}\in\mathbb{R}^{8}$ be the vector of all $1$’s.

Then, the system

(2.13) $M(p)\varphi(p)=\begin{bmatrix}1\\\ p\\\ 0_{4}\end{bmatrix},\quad
M(p)\mathrel{\mathop{:}}=\begin{bmatrix}\mathds{1}^{\top}\\\ V\\\ \Delta(p)\\\
d(p)^{\top}\end{bmatrix}$

has a unique, componentwise nonnegative solution $\varphi(p)$.

However, it can be seen that if $p\in\mathrm{bnd}\,H$ then (2.13), while still
nonsingular, can provide a unique unfeasible solution $\varphi(p)$. This issue
can be overcome by suitably placing zeros on all the columns of $\Delta(p)$ in
(2.12) relative to the indices of the face that $p$ belongs to. In order to do
so, letting $\mathcal{F}$ be the planar face with vertices
$v_{i},v_{j},v_{k},v_{h}$, if it holds that

$(p-v_{i})\times(v_{j}-v_{i})=(v_{k}-v_{i})\times(v_{j}-v_{i}),$

then $p\in\mathcal{F}$ and thus we set

$\Delta(:,[i,j,k,h])=0_{3\times 4}.$

It is a consequence from the 2D case that the corresponding $M(p)$ is
invertible and that the unique solution $\varphi(p)$ to (2.13) is
componentwise nonnegative, such that $\varphi_{r}(p)=0$, for $r=i,j,k,h$.

It is fundamental to point out that (2.11) is a crucial assumption to obtain
feasibility of the unique solution to (2.13). A simple argument is needed to
show that (2.11) is always satisfied, for each $p\in H$, if $H$ is a regular
hexahedron (for example, if $H$ is a cube). However, when $H$ is a general
convex hexahedron, one can easily see that (2.11) may fail. In this case, a
simple change of coordinates can fix this issue. More precisely, if $H$ is
convex, one can always find three linearly independent normals, providing
three planes, such that each vertex in $V$ lies in a different orthant, as
shown below.

######  Lemma 2.5.

Let $H$ be a convex hexahedron. Then, for each $p\in H$ there exists a unit
reference system $R(p)\mathrel{\mathop{:}}=\\{r_{1}(p),r_{2}(p),r_{3}(p)\\}$
such that the sign pattern of $V-p$, in $R(p)$ is given by (2.11).

###### Proof.

Let $F_{i},G_{i},H_{i},i=1,2$ be the three couples of opposite faces of $H$.
If $F_{1},F_{2}$ are parallel, then let $r_{1}$ be any vector parallel to
$F_{1}$. Otherwise, if $F_{1},F_{2}$ are not parallel, let
$l\mathrel{\mathop{:}}=S_{F_{1}}\cap S_{F_{2}}$ be the portion of line at the
intersection of the supporting planes of $F_{1},F_{2}$ and within the
polyhedron defined by the supporting planes of all the other faces. Let $A\in
l$ be arbitrarily chosen. Then, let $r_{1}$ be the unit vector relative to
$p-A$. The remaining vectors $r_{2},r_{3}$ are constructed analogously as
above, according to whether $G_{i},H_{i},i=1,2$ are parallel or not. The set
$\\{r_{1},r_{2},r_{3}\\}$ defined this way is then a basis for
$\mathbb{R}^{3}$ such that $\mathrm{sgn}(V-p)$ is as in (2.11). ∎

In addition, we observe that if $H$ is a nonconvex hexahedron, the conclusion
of Lemma 2.5 may fail depending on the point $p\in H$.

##  3\. Numerical results

We present analytical expressions of moment coordinates over convex and
nonconvex quadrilaterals, and present plots of these coordinates over
quadrilaterals and hexahedra.

###  3.1. Quadrilateral

On solving (2.2) over the biunit square, $Q=\square$, we obtain the analytical
solution for the moment coordinates (identical to mean value coordinates):

$\varphi=\begin{bmatrix}-\frac{x\sqrt{x^{2}-2x+y^{2}-2y+2}-\sqrt{x^{2}+2x+y^{2}-2y+2}-\sqrt{x^{2}-2x+y^{2}+2y+2}+x\sqrt{x^{2}+2x+y^{2}-2y+2}+y\sqrt{x^{2}-2x+y^{2}-2y+2}+y\sqrt{x^{2}-2x+y^{2}+2y+2}}{2\left(\sqrt{x^{2}-2x+y^{2}-2y+2}+\sqrt{x^{2}-2x+y^{2}+2y+2}+\sqrt{x^{2}+2x+y^{2}-2y+2}+\sqrt{x^{2}+2x+y^{2}+2y+2}\right)}\\\
\frac{\sqrt{x^{2}-2x+y^{2}-2y+2}+\sqrt{x^{2}+2x+y^{2}+2y+2}+x\sqrt{x^{2}-2x+y^{2}-2y+2}+x\sqrt{x^{2}+2x+y^{2}-2y+2}-y\sqrt{x^{2}+2x+y^{2}-2y+2}-y\sqrt{x^{2}+2x+y^{2}+2y+2}}{2\left(\sqrt{x^{2}-2x+y^{2}-2y+2}+\sqrt{x^{2}-2x+y^{2}+2y+2}+\sqrt{x^{2}+2x+y^{2}-2y+2}+\sqrt{x^{2}+2x+y^{2}+2y+2}\right)}\\\
\frac{\sqrt{x^{2}-2\,x+y^{2}+2y+2}+\sqrt{x^{2}+2x+y^{2}-2\,y+2}+x\sqrt{x^{2}-2x+y^{2}+2y+2}+x\sqrt{x^{2}+2x+y^{2}+2y+2}+y\sqrt{x^{2}+2x+y^{2}-2y+2}+y\sqrt{x^{2}+2x+y^{2}+2y+2}}{2\left(\sqrt{x^{2}-2x+y^{2}-2y+2}+\sqrt{x^{2}-2x+y^{2}+2y+2}+\sqrt{x^{2}+2x+y^{2}-2y+2}+\sqrt{x^{2}+2x+y^{2}+2y+2}\right)}\\\
\frac{\sqrt{x^{2}-2x+y^{2}-2y+2}+\sqrt{x^{2}+2x+y^{2}+2y+2}-x\sqrt{x^{2}-2x+y^{2}+2y+2}-x\sqrt{x^{2}+2x+y^{2}+2y+2}+y\sqrt{x^{2}-2x+y^{2}-2y+2}+y\sqrt{x^{2}-2x+y^{2}+2y+2}}{2\left(\sqrt{x^{2}-2x+y^{2}-2y+2}+\sqrt{x^{2}-2x+y^{2}+2y+2}+\sqrt{x^{2}+2x+y^{2}-2y+2}+\sqrt{x^{2}+2x+y^{2}+2y+2}\right)}\end{bmatrix}$

As is well known, mean value coordinates are distinct from bilinear finite
element shape functions on the biunit square. However, Wachspress coordinates
are identical to bilinear finite element shape functions which is realized on
solving (2.6). We now present these coordinates on two simple (convex and
nonconvex) quadrilaterals.

######  Example 3.1.

Let us consider the convex quadrilateral $Q$ with vertices

$v_{1}=\begin{bmatrix}0\\\ 0\end{bmatrix},\ v_{2}=\begin{bmatrix}1\\\
0\end{bmatrix},\ v_{3}=\begin{bmatrix}\frac{1}{2}\\\ 4\end{bmatrix},\
v_{4}=\begin{bmatrix}0\\\ 2\end{bmatrix}.$

Figure 4 shows both mean value and Wachspress coordinates for each vertex.
Both coordinates are nonnegative, piecewise affine on the boundary and satisfy
the Kronecker-delta property at the vertices. Furthermore, we obtain the
analytical solution for Wachspress coordinates:

$\varphi=\left[\begin{array}[]{c}\frac{-32\,x^{2}+4\,x\,y+16\,x+y^{2}-10\,y+16}{2\,\left(16\,x-y+8\right)}\\\
\frac{4\,x\,\left(4\,x-y+2\right)}{16\,x-y+8}\\\ \frac{6\,x\,y}{16\,x-y+8}\\\
\frac{y\,\left(8-8\,x-y\right)}{2(16\,x-y+8)}\end{array}\right],$

which are rational functions, with numerator of degree $2$ and denominator of
degree $1$ [4].

(a)

(b)

Figure 4. Plots of (a) moment and (b) Wachspress coordinates on $Q$ for
Example 3.1.

(a)

(b)

Figure 6. Plots of (a) $x$-derivative and (b) $y$-derivative of moment
coordinates on $Q$ for Example 3.1.

(a)

(b)

Figure 8. Plots of (a) $x$-derivative and (b) $y$-derivative of Wachspress
coordinates on $Q$ for Example 3.1.

######  Example 3.2.

Let us consider the nonconvex quadrilateral $Q$ with vertices

$v_{1}=\begin{bmatrix}0\\\ 0\end{bmatrix},\ v_{2}=\begin{bmatrix}2\\\
0\end{bmatrix},\ v_{3}=\begin{bmatrix}1\\\ 4\end{bmatrix},\
v_{4}=\begin{bmatrix}0\\\ 2\end{bmatrix}.$

On a nonconvex quadrilateral, Wachspress coordinates are not valid but moment
(mean value) coordinates are permissible. In Figure 9, moment coordinates are
presented for each vertex and once again we observe that these coordinates
meet the desired properties of generalized barycentric coordinates given in
(1.1). The analytical solutions for moment coordinates in this case are:

$\varphi=\begin{bmatrix}-\frac{2\,\sqrt{x^{2}-4\,x+y^{2}+4}-8\,\sqrt{x^{2}-2\,x+y^{2}-4\,y+5}-4\,\sqrt{x^{2}-2\,x+y^{2}-8\,y+17}-2\,x\,\sqrt{x^{2}-4\,x+y^{2}+4}+4\,x\,\sqrt{x^{2}-2\,x+y^{2}-4\,y+5}+2\,x\,\sqrt{x^{2}-2\,x+y^{2}-8\,y+17}+y\,\sqrt{x^{2}-2\,x+y^{2}-4\,y+5}+y\,\sqrt{x^{2}-2\,x+y^{2}-8\,y+17}}{2\,\left(4\,\sqrt{x^{2}-2\,x+y^{2}-4\,y+5}+2\,\sqrt{x^{2}-2\,x+y^{2}-8\,y+17}-\sqrt{x^{2}-4\,x+y^{2}+4}+\sqrt{x^{2}+y^{2}}\right)}\\\
\frac{2\,x\,\sqrt{x^{2}+y^{2}}-2\,\sqrt{x^{2}+y^{2}}+4\,x\,\sqrt{x^{2}-2\,x+y^{2}-4\,y+5}+2\,x\,\sqrt{x^{2}-2\,x+y^{2}-8\,y+17}-y\,\sqrt{x^{2}-2\,x+y^{2}-4\,y+5}-y\,\sqrt{x^{2}-2\,x+y^{2}-8\,y+17}}{2\,\left(4\,\sqrt{x^{2}-2\,x+y^{2}-4\,y+5}+2\,\sqrt{x^{2}-2\,x+y^{2}-8\,y+17}-\sqrt{x^{2}-4\,x+y^{2}+4}+\sqrt{x^{2}+y^{2}}\right)}\\\
\frac{2\,x\,\sqrt{x^{2}-4\,x+y^{2}+4}-y\,\sqrt{x^{2}-4\,x+y^{2}+4}+2\,x\,\sqrt{x^{2}+y^{2}}+y\,\sqrt{x^{2}+y^{2}}-4\,\sqrt{x^{2}+y^{2}}+2\,y\,\sqrt{x^{2}-2\,x+y^{2}-4\,y+5}}{2\,\left(4\,\sqrt{x^{2}-2\,x+y^{2}-4\,y+5}+2\,\sqrt{x^{2}-2\,x+y^{2}-8\,y+17}-\sqrt{x^{2}-4\,x+y^{2}+4}+\sqrt{x^{2}+y^{2}}\right)}\\\
-\frac{4\,x\,\sqrt{x^{2}-4\,x+y^{2}+4}-y\,\sqrt{x^{2}-4\,x+y^{2}+4}+4\,x\,\sqrt{x^{2}+y^{2}}+y\,\sqrt{x^{2}+y^{2}}-8\,\sqrt{x^{2}+y^{2}}-2\,y\,\sqrt{x^{2}-2\,x+y^{2}-8\,y+17}}{2\,\left(4\,\sqrt{x^{2}-2\,x+y^{2}-4\,y+5}+2\,\sqrt{x^{2}-2\,x+y^{2}-8\,y+17}-\sqrt{x^{2}-4\,x+y^{2}+4}+\sqrt{x^{2}+y^{2}}\right)}\end{bmatrix}.$

Figure 9. Moment coordinates on $Q$ for Example 3.2.

(a)

(b)

Figure 11. Plots of (a) $x$-derivative and (b) $y$-derivative of moment
coordinates on $Q$ for Example 3.2.

###  3.2. Hexahedron

######  Example 3.3.

Let us consider the convex hexahedron $H_{c}$ with vertices

$v_{1}=\begin{bmatrix}1\\\ 2\\\ 1\end{bmatrix},\ v_{2}=\begin{bmatrix}1\\\
2\\\ -1\end{bmatrix},\ v_{3}=\begin{bmatrix}1\\\ 0\\\ -1\end{bmatrix},\
v_{4}=\begin{bmatrix}1\\\ 0\\\ 1\end{bmatrix},v_{5}=\begin{bmatrix}-1\\\ 1\\\
1\end{bmatrix},\ v_{6}=\begin{bmatrix}-1\\\ 1\\\ -1\end{bmatrix},\
v_{7}=\begin{bmatrix}-1\\\ -1\\\ -1\end{bmatrix},\ v_{8}=\begin{bmatrix}-1\\\
-1\\\ 1\end{bmatrix}$

Moment coordinates are nonnegative on the whole hexahedron, piecewise affine
on the boundary edges and satisfy the Kronecker-delta property at the
vertices; moreover, they satisfy the facet-reduction property on each face, as
shown in Figure 12.

Figure 12. Moment coordinates on $H_{c}$ for Example 3.3.

##  4\. Conclusions

In this paper, we introduced moment coordinates, $\varphi(p)$, a new
generalized barycentric coordinates over common finite element geometries.
Generalized barycentric coordinates must satisfy the linear reproducing
conditions and it is desirable that they be nonnegative. We supplemented the
reproducing conditions by moment regularizing equality constraints so that the
linear system of equations yielded $\varphi(p)$ that were unique and
componentwise nonnegative. This idea of adding moment equality conditions goes
back to work on sliding mode control in piecewise smooth dynamical systems
[16, 17], in which one searches for Filippov solutions in presence of dynamics
along discontinuity manifolds of codimension 2 or higher. We showed that by
judicious choice of the regularizing constraints, we could recover mean value
coordinates [6] on nonconvex quadrilaterals as well as Wachspress coordinates
[4] on convex quadrilaterals. Moment coordinates on general convex hexahedron
were constructed by using signed partial distances [17] in tandem with a local
reference system, which ensured that $V-p$ ($V$ is the matrix of vertex
coordinates) had the same sign pattern for all points in the hexahedron.
Theoretical proofs were presented which were supported by analytical
expressions for moment coordinates and their plots over quadrilaterals and
hexahedra.

Extension of the present work to polygonal finite elements in two dimensions
and to other finite element geometries (for example, prisms and pyramids) in
three dimensions are of interest. For a convex polytope with $n$ vertices,
maximum-entropy regularization [22] is a suitable approach to obtain
$\varphi(p)$, but it is not aligned with our objective of obtaining a closed-
form expression for $\varphi(p)$. In keeping with the approach presented in
this paper, it is desirable to devise $n$ linear constraints that yield a
feasible nonnegative solution for $\varphi(p)$. To this end, as a first step
in future work we plan to study the extension of moment coordinates to simple
polygons.

## Acknowledgments

FVD has been supported by REFIN Project, grant number 812E4967. He also
gratefully acknowledges INdAM Research Group GNCS.

## Conflict of interest

The authors declare no conflict of interest.

## References

  * 1  M. S. Floater, Generalized barycentric corodinates and applications, _Acta Numerica_ , 24 (2015), 161—-214
  * 2  D. Anisimov, Barycentric coordinates and their properties, In: K. Hormann and N. Sukumar, _Generalized Barycentric Coordinate in Computer Graphics and Computational Mechanics_ , Boca Raton: Taylor & Francis, CRC Press, 2017, 3—-22
  * 3  K. Hormann and N. Sukumar, _Generalized Barycentric Coordinates in Computer Graphics and Computational Mechanics_ , Boca Raton: Taylor & Francis, CRC Press, 2017
  * 4  E. L. Wachspress, _Rational Bases and Generalized Barycentric_ , Cham: Springer, 2016
  * 5  J. Warren, _Barycentric coordinates for convex polytopes_. _Adv. Comput. Math._ , 6:1 (1996), 97—-108
  * 6  M. S. Floater, Mean value coordinates, _Comput. Aided Geom. Des._ , 20:1 (2003), 19—-27
  * 7  K. Hormann and M. S. Floater, Mean value coordinates for arbitrary planar polygons, _ACM Trans. Graph._ , 25:4 (2006), 1424—-1441
  * 8  M. Arroyo and M. Ortiz, Local maximum-entropy approximation schemes: A seamless bridge between finite elements and meshfree methods, _Int. J. Numer. Method Eng._ , 65:13 (2006), 2167—-2202
  * 9  P. Joshi, M. Meyer, T. DeRose, B. Green and T. Sanocki, Harmonic coordinates for character articulation, _ACM Trans. Graph._ , 26:3 (2007), Article 71, 9 pages
  * 10  N. Sukumar and R. W. Wright, Overview and construction of meshfree basis functions: from moving least squares to entropy approximants, _Int. J. Numer. Method Eng._ , 70:2 (2007), 181—-205
  * 11  K. Hormann and N. Sukumar, Maximum entropy coordinates for arbitrary polytopes, _Comput. Graph. Forum_ , 27:5 (2008), 1513—-1520
  * 12  J. Zhang, B. Deng, Z. Liu, G. Patanè, S. Bouaziz, K. Hormann and L. Liu, _ACM Trans. Graph._ , 33:6 (2014), Article 188, 12 pages
  * 13  J. Tao, B. Deng and J. Zhang, A fast numerical solver for local barycentric coordinates. _Comput. Aided Geom. Des._ , 70 (2019), 46—-58
  * 14  D. Anisimov, D. Panozzo and K. Hormann, Blended barycentric coordinates. _Comput. Aided Geom. Des._ , 52–53 (2017), 205—-216
  * 15  C. Deng, Q. Chang and K. Hormann, _Comput. Aided Geom. Des._ , 79 (2020), Article 101861
  * 16  L. Dieci and F. Difonzo, The moments sliding vector field on the intersection of two manifolds, _J. Dyn. Differ. Equations_. 29:1 (2017), 169—-201
  * 17  L. Dieci and F. Difonzo, On the inverse of some sign matrices and on the moments sliding vector field on the intersection of several manifolds: nodally attractive case, _J. Dyn. Differ. Equations_. 29:4 (2017), 1355—-1381.
  * 18  V. Klee, R. Ladner, R. Manber, Signsolvability revisited, _Linear Algebra Appl._ , 59 (1984), 131—-157
  * 19  R. A. Brualdi, B. L. Shader, B. Bollobas, _Matrices of Sign-Solvable Linear Systems_ , Cambridge Tracts in Mathematics, Cambridge University Press, 1995
  * 20  D. P. Flanagan and T. Belytschko, A uniform strain hexahedron and quadrilateral with orthogonal hourglass control, _Int. J. Numer. Method Eng._ , 17:5 (1981), 679—-706
  * 21  M. Floater, A. Gillette, N. Sukumar, Gradient bounds for Wachspress coordinates on polytopes, _SIAM J. Numer. Anal._ , 52:1 (2014), 515—-532
  * 22  N. Sukumar, Construction of polygonal interpolants: a maximum entropy approach, _Int. J. Numer. Method Eng._ , 61:12 (2004), 2004—-2181