Mapping an arbitrary 4 node quadrilateral element in space (x,y,z) to unit square

Question

I need to integrate a function G(x,y,z) over an arbitrary 4 node quadrilateral element (4 nodes xi,yi,zi on a same plane). I know how to handle it when the 4 nodes (xi,yi) are on the X-Y plane, i.e., the isoparametric mapping. But when the nodes' coordinates are 3 dimensional (xi,yi,zi), is there a direct way like a shape function or jacobian matrix to do that?

Answer 1

To map the unit square $(u, v)$, $0 \le u, v \le 1$ to an arbitrary quadrilateral $(x_i, y_i, z_i)$, $i = 0, 1, 2, 3$, you can use bilinear interpolation: $$\begin{aligned} x(u,v) &= (1 - u) (1 - v) x_0 + u (1 - v) x_1 + (1 - u) v x_3 + u v x_2 \\ y(u,v) &= (1 - u) (1 - v) y_0 + u (1 - v) y_1 + (1 - u) v y_3 + u v y_2 \\ z(u,v) &= (1 - u) (1 - v) z_0 + u (1 - v) z_1 + (1 - u) v z_3 + u v z_2 \\ \end{aligned}$$ or, equivalently, in expanded form $$\begin{aligned} x(u,v) &= x_0 + (x_1 - x_0) u + (x_3 - x_0) v + (x_0 - x_1 + x_2 - x_3) u v \\ y(u,v) &= y_0 + (y_1 - y_0) u + (y_3 - y_0) v + (y_0 - y_1 + y_2 - y_3) u v \\ z(u,v) &= z_0 + (z_1 - z_0) u + (z_3 - z_0) v + (z_0 - z_1 + z_2 - z_3) u v \\ \end{aligned}$$ The latter form is useful, if you just choose to do the integration in 3D, as you can directly see the limits and directional derivatives along the $u$ and $v$ axes. It should be a relatively simple change in variables.

Bilinear interpolation is valid at least when the quadrilateral is planar and non-self-intersecting. If it is planar and self-intersecting (bowtie, ⧓), then depending on the labeling of the four vertices, the interpolation may be incorrect/outside the actual figure. If the quadrilateral is nonplanar, the interpolation chooses a specific curvature of the surface that only depends on the 3D vertices. That may or may not be desirable.

The inverse exists unless the quadrilateral is degenerate somehow. Numerically, you pick the two axes ($x,y$, $x,z$, or $y,z$) maximizing the projected area of the quadrilateral to minimize the rounding errors, and then solve the pair of equations above for $u$ and $v$.

---

Computationally, planar surfaces are typically split into triangles instead of quadrilaterals, using barycentric coordinates ($(u, v, w)$, $0 \le u, v, w \le 1$, $ u + v + w = 1$). The center of the triangle is at $(1/3, 1/3, 1/3)$. Note that since $u$ and $v$ completely define the third barycentric coordinate $w$, $w$ is usually omitted. If we number the three vertices so that $0: (0, 0, 1)$, $1: (1, 0, 0)$, and $2: (0, 1, 0)$, in barycentric coordinates $(u, v, w)$, then $$\left\lbrace \begin{aligned} x(u, v, w) &= x_0 w + x_1 u + x_2 v \\ y(u, v, w) &= y_0 w + y_1 u + y_2 v \\ z(u, v, w) &= z_0 w + z_1 u + z_2 v \\ \end{aligned}\right. , \quad \left\lbrace\begin{aligned} 0 \le u &\le 1 \\ 0 \le v &\le 1 \\ 0 \le w &\le 1 \\ u + v + w &= 1 \\ \end{aligned}\right.$$ or, equivalently, omitting $w$ coordinate, $$\left\lbrace \begin{aligned} x(u, v) &= x_0 + (x_1 - x_0) u + (x_2 - x_0) v \\ y(u, v) &= y_0 + (y_1 - y_0) u + (y_2 - y_0) v \\ z(u, v) &= z_0 + (z_1 - z_0) u + (z_2 - z_0) v \\ \end{aligned} \right. , \quad \left\lbrace \begin{aligned} 0 \le u \le 1 - v \le 1 \\ 0 \le v \le 1 - u \le 1 \\ \end{aligned} \right.$$ Because triangles are always planar in 3D, if the triangle is not degenerate (is not a line segment or a point; has nonzero area), this does have an inverse. Again, pick two coordinate axes, and solve for $u$ and $v$. For example, picking $x$ and $y$, the inverse solution is $$\left\lbrace \begin{aligned} u &= \displaystyle \frac{ (y_2 - y_0) x + (x_0 - x_2) y + x_2 y_0 - x_0 y_2 }{ (x_1 - x_0)(y_2 - y_0) - (x_2 - x_0)(y_2 - y_0) } \\ v &= \displaystyle \frac{ (y_0 - y_1) x + (x_1 - x_0) y + x_0 y_1 - x_1 y_0 }{ (x_1 - x_0)(y_2 - y_0) - (x_2 - x_0)(y_2 - y_0) } \\ \end{aligned} \right.$$ noting that the divisor is twice the area of the triangle.