Say we have 4 non-planar points in a 3D space (but preferably it should work in any dimension, so no cross products), $a,b,c,d$, and these points define a non-flat area. A point on that area may be computed with $$p_{s,t}=(1-s)(1-t)a+(1-s)tb+s(1-t)c+std$$
How is the surface area of such a figure computed?
I have very little justification for this, but a voice in my head is telling me that this might work: $$\int_{0}^{1}\int_{0}^{1}\frac{1}{d\cdot n_{s,t}}dsdt$$ where $d=unit(\left ( \left ( b-a \right )\times \left ( d-c \right ) \right )\times \left ( \left ( c-a \right )\times \left ( d-b \right ) \right ))$ and $n_{s,t}$ is the normal at point $p_{s,t}$
EDIT: A closed form solution, if possible, would be awesome.
In $3$d the area of a surface defined parametrically by $\mathbf r(s,t): D\subset \Bbb R^2 \to \Bbb R^3$ is given by $$S = \int_D \|\partial_s \mathbf r \times \partial_t \mathbf r\|\ dA = \int_{t_0}^{t_1}\int_{s_0}^{s_1} \|\partial_s \mathbf r \times \partial_t \mathbf r\|\ dsdt$$ where $D = \{(s,t)\mid s_0 \le s\le s_1 \text{ and } t_0 \le t\le t_1\}$.
What you don't learn in most multivariable courses is that the generalization of this formula for surfaces defined in $\Bbb R^n$ for any $n\ge 2$ is $$S = \int_D \|\partial_s \mathbf r \wedge \partial_t \mathbf r\|\ dA$$