Calculate surface normal and area for a non-planar quadrilateral Given the four coordinates of the vertices, what is the best possible approximation to calculate surface area and outward normal for a quad? 
I currently join the midpoints of the sides, thus dividing the quad into four triangles and a parallelogram, and find the normal to the plane of the parallelogram, but that doesn't give me correct results when the original four vertices aren't coplanar.
Is it possible to use isoparametric transformations for this purpose?  (This calculation is a part of a bigger FEM code) 
 A: Lets examine the quadrilateral as a linearly interpolated surface. If the vertices are
$$\bbox{
\vec{p}_1 = \left [ \begin{matrix} x_1 \\ y_1 \\ z_1 \end{matrix} \right ]
} , \quad \bbox {
\vec{p}_2 = \left [ \begin{matrix} x_2 \\ y_2 \\ z_2 \end{matrix} \right ]
} , \quad \bbox {
\vec{p}_3 = \left [ \begin{matrix} x_3 \\ y_3 \\ z_3 \end{matrix} \right ]
} , \quad \bbox {
\vec{p}_4 = \left [ \begin{matrix} x_4 \\ y_4 \\ z_4 \end{matrix} \right ]
}$$
and we parametrise the surface using a vector-valued function $\vec{p}(u, v)$, with $0 \le u \le 1$, $0 \le v \le 1$, then
$$\bbox{ \vec{p}(u, v) = (1 - v)\biggr( (1-u) \vec{p}_1 + u \vec{p}_2 \biggr) + v \biggr( (1-u) \vec{p}_3 + u \vec{p}_4 \biggr) }$$
i.e.
$$\bbox{ \vec{p}(u, v) = \vec{p}_1 + u \left( \vec{p}_2 - \vec{p}_1 \right) + v \left( \vec{p}_3 - \vec{p}_1 \right) + u v \left( \vec{p}_4 - \vec{p}_3 + \vec{p}_1 - \vec{p}_2 \right) }$$
If the four vertices are nonplanar, the surface is curved. In all cases, the surface passes through 
$$\bbox{ \vec{p}(0.5, 0.5) = \frac{ \vec{p}_1 + \vec{p}_2 + \vec{p}_3 + \vec{p}_4 }{4} }$$
The surface tangent vectors are
$$\bbox{ \begin{aligned}
\vec{p}_u (u , v) &= \frac{\partial \vec{p}(u, v)}{\partial u} = \vec{p}_2 - \vec{p}_1 + v \left( \vec{p}_4 - \vec{p}_3 - \vec{p}_2 + \vec{p}_1 \right)
\\
\vec{p}_v (u, v) &= \frac{\partial \vec{p}(u, v)}{\partial v} = \vec{p}_3 - \vec{p}_1 + u \left( \vec{p}_4 - \vec{p}_3 - \vec{p}_2 + \vec{p}_1 \right) \\
\end{aligned} }$$
and the exact area of this surface is
$$\bbox{ A = \int_0^1 d u \int_0^1 d v \; \vec{p}_u (u, v) \cdot \vec{p}_v (u, v) = \frac{ \left( \vec{p}_1 + \vec{p}_2 - \vec{p}_3 - \vec{p}_4 \right)\cdot\left( \vec{p}_1 - \vec{p}_2 + \vec{p}_3 - \vec{p}_4 \right) }{4} }$$
The surface normal vector is $\vec{n}(u, v)$,
$$\bbox{ \vec{n}(u , v) = \frac{\partial \vec{p}(u, v)}{\partial u} \times \frac{\partial \vec{p}(u, v)}{\partial v} }$$
and the mean (expected value) of the normal vector is
$$\bbox{ \langle\vec{n}\rangle = \int_0^1 d u \int_0^1 d v \; \vec{n}(u ,v) = \left [ \begin{matrix}
\frac{(y_4 - y_1)(z_3 - z_2) - (z_4 - z_1)(y_3 - y_2)}{2} \\
\frac{(z_4 - z_1)(x_3 - x_2) - (x_4 - x_1)(z_3 - z_2)}{2} \\
\frac{(x_4 - x_1)(y_3 - y_2) - (y_4 - y_1)(x_3 - x_2)}{2} \\
\end{matrix} \right] = \vec{n}\left(\frac{1}{2} , \frac{1}{2}\right) }$$ 
If we consider a flow through the surface, we could assume the flow is perpendicular to the surface normal. Then, the effective surface is the area of the surface when projected to a plane perpendicular to the unit normal vector $\hat{n}$ used,
for example
$$\hat{n} = \frac{\langle\vec{n}\rangle}{\left\lVert\langle\vec{n}\rangle\right\rVert}$$
$$A_\perp = \int_0^1 d u \int_0^1 d v \; \Biggr( \vec{p}_u (u, v) - \hat{n} \biggr(\hat{n}\cdot\vec{p}_u (u,v) \biggr) \Biggr) . \Biggr( \vec{p}_v (u, v) - \hat{n} \biggr(\hat{n}\cdot\vec{p}_v (u,v) \biggr) \Biggr)$$
It turns out that for this choice of $\hat{n}$, $A_\perp = A$.
Therefore, choosing
$$\bbox[#ffffef, 1em]{ \vec{n} = \left [ \begin{matrix}
(y_4 - y_1)(z_3 - z_2) - (z_4 - z_1)(y_3 - y_2) \\
(z_4 - z_1)(x_3 - x_2) - (x_4 - x_1)(z_3 - z_2) \\
(x_4 - x_1)(y_3 - y_2) - (y_4 - y_1)(x_3 - x_2) \end{matrix} \right ] = \; \bigr( \vec{p}_4 - \vec{p}_1 \bigr) \times \bigr( \vec{p}_3 - \vec{p}_2 \bigr) }$$
$$\bbox{\hat{n} = \frac{\vec{n}}{\left\lVert\vec{n}\right\rVert} = \frac{\vec{n}}{\sqrt{\vec{n} \cdot \vec{n}}}}$$
$$\bbox[#ffffef, 1em]{ A = \frac{ \left( \vec{p}_1 + \vec{p}_2 - \vec{p}_3 - \vec{p}_4 \right)\cdot\left( \vec{p}_1 - \vec{p}_2 + \vec{p}_3 - \vec{p}_4 \right) }{4} }$$
is an obvious choice, in my opinion.
A: The area formula in the posted answer is wrong.  Putting in a square, for example, produces zero area because its sides are normal.   
Using the same framework as Nominal Animal's answer, I believe the correct area formula is $$A=\int_0^1 \int_0^1 |p_u(u,v)\times p_v(u,v)|\,du\,dv,$$ which is a much more complex integral, that I have been unable to reduce to a simple formula.  But at least it worked on some coplanar test cases.  
A: Adrian is right, and Nominal Animals answer has the incorrect surface integral. The correct integral is
$$\begin{aligned}
A & = \int_0^1 \int_0^1 \left\lVert \vec{p}_u \times \vec{p}_v \right \rVert \, d u \, d v \\
~ & = \int_0^1 \int_0^1 \sqrt{ u^2 \, C_{uu} + 2 \, u \, v \, C_{uv} + v^2 \, C_{vv} + 2 \, u \, C_{u} + 2 \, v \, C_{v} + C_{0} } \, d u \, d v \\
\end{aligned}$$
where numbering the vertices such that $1$ and $4$ are diagonally opposite (as are $2$ and $3$),
$$\begin{aligned}
C_{uu} & = ( Y_\Delta Z_2 - Z_\Delta Y_2 )^2
         + ( Z_\Delta X_2 - X_\Delta Z_2 )^2
         + ( X_\Delta Y_2 - Y_\Delta X_2 )^2 \\
C_{vv} & = ( Y_\Delta Z_3 - Z_\Delta Y_3 )^2
         + ( Z_\Delta X_3 - X_\Delta Z_3 )^2
         + ( X_\Delta Y_3 - Y_\Delta X_3 )^2 \\
C_{uv} & = ( Y_\Delta Z_3 - Z_\Delta Y_3 )( Z_\Delta Y_2 - Y_\Delta Z_2 ) \\
~      & + ( Z_\Delta X_3 - X_\Delta Z_3 )( X_\Delta Z_2 - Z_\Delta X_2 ) \\
~      & + ( X_\Delta Y_3 - Y_\Delta X_3 )( Y_\Delta X_2 - X_\Delta Y_2 ) \\
C_{u}  & = ( Y_2 Z_\Delta - Z_2 Y_\Delta )( Y_2 Z_3 - Z_2 Y_3 ) \\
~      & + ( Z_2 X_\Delta - X_2 Z_\Delta )( Z_2 X_3 - X_2 Z_3 ) \\
~      & + ( X_2 Y_\Delta - Y_2 X_\Delta )( X_2 Y_3 - Y_2 X_3 ) \\
C_{v}  & = ( Y_\Delta Z_3 - Z_\Delta Y_3 )( Y_2 Z_3 - Z_2 Y_3 ) \\
~      & + ( Z_\Delta X_3 - X_\Delta Z_3 )( Z_2 X_3 - X_2 Z_3 ) \\
~      & + ( X_\Delta Y_3 - Y_\Delta X_3 )( X_2 Y_3 - Y_2 X_3 ) \\
C_{0}  & = ( Y_2 Z_3 - Z_2 Y_3 )^2 
         + ( Z_2 X_3 - X_2 Z_3 )^2 
         + ( X_2 Y_3 - Y_2 X_3 )^2  \\
\end{aligned}$$
where $C_{uu} \ge 0$, $C_{vv} \ge 0$, and $C_{0} \ge 0$; and
$$\begin{aligned}
X_\Delta & = x_4 - x_3 - x_2 + x_1 = (x_4 - x_3) - (x_2 - x_1) \\
X_2      & = x_2 - x_1 \\
X_3      & = x_3 - x_1 \\
Y_\Delta & = y_4 - y_3 - y_2 + y_1 = (y_4 - y_3) - (y_2 - y_1) \\
Y_2      & = y_2 - y_1 \\
Y_3      & = y_3 - y_1 \\
Z_\Delta & = z_4 - z_3 - z_2 + z_1 = (z_4 - z_3) - (z_2 - z_1) \\
Z_2      & = z_2 - z_1 \\
Z_3      & = z_3 - z_1 \\
\end{aligned}$$
but I am not sure if there is a feasible algebraic expression for the integral.  Numerically it is a simple matter, as the expression in the square root is numerically stable and just a bivariate quadratic function.

To calculate the cross sectional area in some direction $\hat{n}$, $\lVert\hat{n}\rVert = 1$, for example a flow in a specific direction through the surface, then
$$A_\hat{n} = \int_0^1 \int_0^1 \hat{n} \cdot \left( \vec{p}_u \times \vec{p}_v \right) \, d u \, d v$$
If $\hat{n} = (n_x, n_y, n_z)$, $n_x^2 + n_y^2 + n_z^2 = 1$ and $\hat{n}$ is in the same halfspace as the surface normal (i.e, $\hat{n} \cdot \left( \vec{p}_u \times \vec{p}_v \right) \ge 0$ for $0 \le u \le 1$, $0 \le v \le 1$, then
$$\begin{aligned}
A_\hat{n} & = \frac{n_x}{2} \bigl( (y_4 - y_1) (z_3 - z_2) - (z_4 - z_1) (y_3 - y_2) \bigr) \\
~         & + \frac{n_y}{2} \bigl( (z_4 - z_1) (x_3 - x_2) - (x_4 - x_1) (z_3 - z_2) \bigr) \\
~         & + \frac{n_z}{2} \bigl( (x_4 - x_1) (y_3 - y_2) - (y_4 - y_1) (x_3 - x_2) \bigr) \\
~         & = \frac{1}{2} \hat{n} \cdot \bigl( (\vec{p}_4 - \vec{p}_1) \times (\vec{p}_3 - \vec{p}_2) \bigr) \\
\end{aligned}$$
i.e. half the dot product between the unit flow direction vector and the cross product of the quadrilateral diagonals.  But do remember that this only applies if $\hat{n}$ is and the quadrilateral surface normal are in the same halfspace within the entire surface.
A: This beautiful answer here is just wrong.
$$\bbox{ A = \int_0^1 d u \int_0^1 d v \; \vec{p}_u (u, v) \cdot \vec{p}_v (u, v) = \frac{ \left( \vec{p}_1 + \vec{p}_2 - \vec{p}_3 - \vec{p}_4 \right)\cdot\left( \vec{p}_1 - \vec{p}_2 + \vec{p}_3 - \vec{p}_4 \right) }{4} }$$
When I was faced with the same problem I divided the nonplanar quadrilateral into four triangles, all sharing the midpoint as one of the vertices. The area of a triangle is half the length of the cross product of any two of the sides, so it's easy to calculate. And the estimate you get is about $12\% \pm 10\%$ too large for randomly generated quadrilaterals. The more flat they are the better this estimate becomes.

