Nonnegative moment coordinates on finite element geometries

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:

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

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:} = {\rm{conv}}\, \{v_i\}$, is the set of all the nonnegative solutions $\varphi(p)$ to the full-rank, underdetermined system

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}{{\mathbb{1}} ^\top \\ V}\end{bmatrix}$. It is straightforward to see that $g$ has a specific sign pattern, which can be computed by

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

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

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 most natural 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\overline{v_1v_2}$, so that we expect

to be the unique solution to

But since $\alpha = \frac{\|p-v_1\|}{\|v_2-v_1\|}$, we get

from which, if $d_1 \mathrel{\mathop:} = \|p-v_1\|, d_2 \mathrel{\mathop:} = -\|p-v_2\|$, the above holds true (see Figure 1a, 1b). Repeating the same argument for the remaining edges, we see that defining

provides (2.2) to be nonsingular and its unique solution to be feasible 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

and

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

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\neq0$ on $\mathbb{R}^3$. Moreover, on the boundary of $[0, 1]^3$ it attains the following values:

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) = \frac12$, 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_1$ out and considering the triangle $v_2, v_3, v_4$, then the corresponding triangular barycentric coordinates $\tau^1(p)$ are computed for each $p\in Q$ as

Applying Cramer's rule to (2.2) for computing $\varphi_1(p)$ provides

where $\mathbf{n} \mathrel{\mathop:} = \left[ \begin{matrix}{A_{1} \\ -A_{2} \\ A_{3} \\ -A_{4}}\end{matrix} \right]$, $A_{i} \mathrel{\mathop:} = \lvert{\mathrm{area}({\rm{conv}}\, \{v_j\, :\, j\neq i\})}\rvert$, and $\langle\mathbf{n}\rangle = \ker\left[ \begin{matrix}{{\mathbb{1}} ^\top \\ V}\end{matrix} \right]$. Letting $\hat\tau^1(p) \mathrel{\mathop:} = \left[ \begin{matrix}{\tau^1_2(p) \\ \tau^1_3(p) \\ \tau^1_4(p)}\end{matrix} \right]\in \mathbb{R}^3$ and observing that $\mathbf{n}^1(p) \mathrel{\mathop:} = A_1\left[ \begin{matrix}{-1 \\ \hat\tau^1(p)}\end{matrix} \right]$ is such that

we have

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

2.2.   Convex quadrilaterals

In this section we prove that on a convex quadrilateral $Q$ (see Figure 1a), akin to mean value 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

where

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

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 $\overline{v_iv_{i+1}}$ and $h_i(p)$ be the distance from $p$ to the edge $\overline{v_iv_{i+1}}$, we define the weights

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

It then follows that

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

and using planar geometry, it is readily shown that

Thus, letting

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

We point out that $p\in Q$, $p\in\overline{v_iv_{i+1}}$ if and only if $h_i(p) = 0$. Therefore, $p\in\overline{v_{i-1}v_i}\cup\overline{v_iv_{i+1}}$ if and only if $\rho_i(p) = 0$.

2.3.   One-dimensional case

Moment conditions 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

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

The system above is easily seen to be invertible, and since $\left[ \begin{matrix}{0\\0}\end{matrix} \right]$ is in the convex hull of the points $\left[ \begin{matrix}{x_1-x\\|x_1-x|}\end{matrix} \right], \left[ \begin{matrix}{x_2-x\\-|x_2-x|}\end{matrix} \right], \left[ \begin{matrix}{x_3-x\\|x_3-x|}\end{matrix} \right]$, then its unique solution $\varphi(x)$ is necessarily feasible.

In the following we extend such an argument to any number $n$ of distinct points $\{x_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

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

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

with $\xi\neq\eta$, then $\left[ \begin{matrix} A \\ d^{\top} \end{matrix} \right]$ has rank $n+1$.

Let $n\geq4$, $V \mathrel{\mathop:} = \left[ \begin{matrix}{x_1 \ \ \ \cdots \ \ \ x_n}\end{matrix} \right]$, with $0\leq x_1 < x_2 < \ldots < x_n\leq1$. Then, given $p\in[0, 1]$, we seek for a feasible solution $\varphi(p)$ of the underdetermined system

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

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

We prove the following.

Proposition 2.3. Let $p\in[0, 1]$. For all $n\geq4$ the system

is invertible and its unique solution is feasible.

Proof. We proceed by induction on $n$. Observe that the base of induction is the case $n = 4$. If so,

Assume $p = x_1$, let $\alpha\neq0$ such that

Therefore

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

Again

and Lemma 2.2 proves the claim.

Let the result hold true for $n-1$, and let us define

Let $p\in(x_1, x_2)$. Then there exist $\alpha, \beta\in(0, 1)$ such that

and

Therefore, assuming $n$ to be even, $\varphi_{12} \mathrel{\mathop:} = \left[ {\begin{array}{*{20}{c}}{1-\alpha }\\{ \alpha }\\{ 0 }\\{ \vdots }\\{ 0}\end{array}} \right], \ \varphi_{1n} \mathrel{\mathop:} = \left[ {\begin{array}{*{20}{c}}{1-\beta }\\{ 0 }\\{ \vdots }\\{ 0 }\\{ \beta}\end{array}} \right]$ are solutions to

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.   Flat-faced 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

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

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

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

Finally, let

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

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

Then, the system

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

However, it can be seen that if $p\in\partial 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

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

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$) that is centered at $[1/2, 1/2]^\top$, we obtain the analytical solution for the moment coordinates (identical to mean value coordinates):

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

Figure 3 shows both mean value and Wachspress coordinates for each vertex, {while their partial derivatives are shown in Figures 4 and 5, respectively}. 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:

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

Example 3.2. Let us consider the nonconvex quadrilateral $Q$ with vertices

On a nonconvex quadrilateral, Wachspress coordinates are not valid but moment (mean value) coordinates are permissible. In Figures 6 and 7, moment coordinates and their partial derivatives, respectively, 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:

Moment coordinates are also well defined and provide the expected solutions in the case of a degenerate quadrilateral that has two consecutive collinear sides, as we show in the next example.

Example 3.3. Let us consider the convex quadrilateral $Q$ with vertices

so that $Q$ has the shape of a triangle.

In Figure 8, moment coordinates are presented for each vertex, whose analytical expressions are given by

As expected from the geometry of the example, the last component $\varphi_4$ is a linear function.

Remark 3.4. It is worth noticing that there is a potential ill-conditioning of the linear system (2.2), but this is a mere reflection of the fact that the barycentric coordinates are not necessarily stable under perturbation. Indeed, this is not a peculiarity of the system formulation.

For example, let us consider the nonconvex quadrilateral $Q_1$ with vertices

Perturbing randomly $Q_1$ to order $10^{-9}$ produces a variation of order $10^{-2}$ in the barycentric coordinates of the point $P = \left[ {\begin{array}{*{20}{c}}{ 0 }\\{ \frac12 10^{-8}} \end{array}} \right]$.

Further, let us consider the rectangle $Q_2$ with vertices

where $\frac{\overline{v_1v_4}}{\overline{v_1v_2}} = 10^{-8}$, so that $\overline{v_1v_4}\ll \overline{v_1v_2}$. With $P = \left[ {\begin{array}{*{20}{c}} {\frac12 }\\{ \frac12 10^{-8}} \end{array}} \right]$, one obtains $\varphi = \frac14{\mathbb{1}}$; however, a random perturbation of the vertices of order $10^{-8}$ creates a variation in the barycentric coordinates of order $1$, regardless of the inexact arithmetic.

3.2.   Hexahedron

Example 3.5. Let us consider the convex hexahedron $H_c$ with vertices

Moment coordinates are nonnegative on the whole hexahedron, piecewise affine on the boundary edges and satisfy the Kronecker-delta property at the vertices; in addition, they satisfy the facet-reduction property on each face, as shown in Figure 9. Furthermore, they provide results that are distinct from 3D mean value coordinates.

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.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

FVD has been supported by REFIN Project, grant number 812E4967, and by INdAM-GNCS 2023 Project, grant number CUP$\_$E53C22001930001.

Conflict of interest

The authors declare no conflicts of interest.