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.



Smooth unit vector field on a tetrahedron to interpolate vertex constraints
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](http://onlinelibrary.wiley.com/doi/10.1002/nme.4366/abstract). For the theory you can look at [Optimal a priori discretization error bounds for geodesic finite elements](http://link.springer.com/article/10.1007%2Fs10208-014-9230-z).


