Smooth unit vector field on a tetrahedron to interpolate vertex constraints

Question

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.

Answer 1

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.