# 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.

• Consider a u for which the norm of the gradient is constant in a small neighbourhod of the vertices and for which u is constant everywhere else. Then the energy is of order eps^3*eps^(-2)=eps, and hence can get arbitrarily small. Usually the b.c. are on the whole boundary and not just at the vertices. The functional usually is considered on H^1 for which point eval.do not make sense. Dec 17, 2014 at 7:45
• Is $$(\nabla v)^2:=\sum_{i=1}^3\sum_{j=1}^3 \left|\frac{\partial v_i}{\partial x_j}\right|^2?$$ Dec 17, 2014 at 8:27
• Both remarks are correct. This means that $E_{\inf}=0$, which indeed is not physically meaningful. One could save smoothness mathematically by using $(\nabla v)^4$. In the application I have in mind, thus the finite element shape functions alone determine the variation of $v$. Thank you for the clarification. Dec 17, 2014 at 11:21
• So if I understand you correctly you want to minimize the energy only roughly but in a way such that you can compute the energy and its derivative by $u_i$? Dec 17, 2014 at 13:01

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.
• Thank you for these useful references. In case of all $u_i$ in one hemisphere, the common approximation schemes work quite well, and I realise that my approach is closely related to Riemannian averaging. The question is mainly about how to get reasonable estimates of $E_{\min}$ in the other case where the $u_i$ are not on a common hemisphere. Dec 16, 2014 at 0:29
• The Riemannian average does in most cases also work when the points do not lie on a common hemisphere. its just that, to my knowledge there is no proof for that. The only case I know where it fails is if $u_i=-u_j$ for some $i,j$ Dec 16, 2014 at 7:05
• Thank you very much for the offer. Of course, I do have some numerical code to find a solution myself. In the 2D case if all products $n_k \cdot u_i>0$ ($n_k$ is the $k$-th outer normal of the triangle and $i,k=1,2,3$) one of the rotations along the triangle sides have to follow the larger angle connecting the two vertex $u_i$'s. There are three possibilities. In the 3D case the situation is similar with six possibilities. The solution to the problem I'm looking for should avoid doing numerical interpolation because in the end I also need the $\partial E_{\min}/\partial u_i$'s. Dec 16, 2014 at 12:55
• By 2D I actually meant a map from a two dimensional domain to $\mathbb{S}^2$. Considering the comment above I now understand the results of my numerical computation. Dec 17, 2014 at 8:32