Smooth unit vector field on a tetrahedron to interpolate vertex constraints For a tetrahedron $T\subset \mathbb{R}^3$ with vertices $r_i\in \mathbb{R}^3$ , $i=1,\ldots,4$, and unit vectors $u_i\in \mathbb{S}^2$ at each vertex $i=1,\ldots,4$ 
consider the (energy) functional
$$
E: C^2(T,\mathbb{S}^2) \to \mathbb{R},~~~ v\mapsto \int \limits _T (\nabla v)^2\, dV~=~ -\int \limits _T v\cdot  \Delta v \, dV
$$
I'm interested whether it is possible to find
an analytical solution of the variational problem 
$$
      E(v) = \min, ~~v(r_i)=u_i,~i=1,\ldots,4.
$$
as a function of the $u_i$.
The unit vector field to such a solution describes a smooth transition between the vertex
constraints, and $v(T)$ is probably the spherical quadrilateral of least surface connecting the $u_i$.
Already an approximation that varies $C^1$-smoothly as a function of the $u_i$ (also for large angles between the $u_i$) and does not underestimate $E_{\rm min}$ would be interesting.
Also a proof that $E_{\min}(u_i)$ is dicontinuous or not differentiable (maybe based on the fact that $v(T)$ has to switch from one to another spherical quadrilateral) would be great.
The problem originates in numerical micromagnetism where $v$  is the local magnetization direction and $E$ the exchange energy in a finite element cell. The commonly used discretization schemes either tend to underestimate $E$ (especially when using the $ \mathbb{R}^3$ instead of $\mathbb{S}^2$ version of $\Delta v$), or have discontinuities for large angles between the $u_i$. I'm currently developing an improved approximation based on the solution for the one-dimensional problem, but would prefer a three-dimensional exact solution if it exists. One idea is to extend the problem to a quarternion field on $T$ for which an analytical solution might exist, although my literature survey so far was unsuccessful.
 A: I dont think there is an explizit expression for the solution. A $C^\infty$ approximation is
$$v(x):=\mathop{\text{argmin}}_{y\in S^2} \sum_{i=1}^4 \lambda_i(x) d^2(y,u_i)$$
where $\lambda_i$ are the barycentric coordinates and $d$ is the spherical distance. This function is well-defined if $(u_i)_{i=1}^4$ lie in the same hemisphere. $v(x)$ is then also called the Riemannian average of $u_i$ with weights $\lambda_i$. You can compute $v(x)$ quite fast by using gradient descent, i.e. the fixpointiteration $y \mapsto \exp_y(\sum_{i=1}^4 \lambda_i(x) \exp^{-1}_y(u_i))$,
where $\exp$ is the exponential map, for the sphere we have the explicit formula
$$\exp_x(y)=\cos(|y|)x+\sin(|y|)\frac{y}{|y|}$$
$$\exp_x^{-1}(y)=\frac{\arccos(x^Ty}{\sqrt{1-(x^Ty)^2}}(y-(x^Ty)x).$$
If you want to solve the optimization problem numerically you can take a look at
Geodesic finite elements on simplicial grids. For the theory you can look at Optimal a priori discretization error bounds for geodesic finite elements.
