[Math] Calculate surface normal and area for a non-planar quadrilateral

areafinite element methodgeometryquadrilateralsurfaces

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)

Best Answer

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.