Computing cross fields -- A PDE approach based on the Ginzburg-Landau theory

TL;DR

This paper introduces a PDE-based method using Ginzburg-Landau theory to compute crossfields on surfaces, aiding in quadrangular mesh generation by identifying well-distributed critical points.

Contribution

It applies Ginzburg-Landau functional to crossfield computation on arbitrary surfaces, providing a novel PDE approach validated through finite element discretization.

Findings

Validated on the unit disk with analytical comparison

Successfully applied to the 2-sphere with critical points at anticube vertices

Demonstrated the method's effectiveness for mesh generation

Abstract

This paper proposes a method to compute crossfields based on the Ginzburg-Landau theory. The Ginzburg-Landau functional has two terms: the Dirichlet energy of the distribution and a term penalizing the mismatch between the fixed and actual norm of the distribution. Directional fields on surfaces are known to have a number of critical points, which are properly identified with the Ginzburg-Landau approach: the asymptotic behavior of Ginzburg-Landau problem provides well-distributed critical points over the 2-manifold, whose indices are as low as possible. The central idea in this paper is to exploit this theoretical background for crossfield computation on arbitrary surfaces. Such crossfields are instrumental in the generation of meshes with quadrangular elements. The relation between the topological properties of quadrangular meshes and crossfields are hence first recalled. It is then shown that a crossfield on a surface can be represented by a complex function of unit norm with a number of critical points, i.e., a nearly everywhere smooth function taking its values in the unit circle of the complex plane. As maximal smoothness of the crossfield is equivalent with minimal energy, the crossfield problem is equivalent to an optimization problem based on Ginzburg-Landau functional. A discretization scheme with Crouzeix-Raviart elements is applied and the correctness of the resulting finite element formulation is validated on the unit disk by comparison with an analytical solution. The method is also applied to the 2-sphere where, surprisingly but rightly, the computed critical points are not located at the vertices of a cube, but at those of an anticube.

Figures

Equations

$$\chi =2-2g-b.$$

$$\chi =n-n_e+n_f$$

$$n_{evf}{n}_f=2(n_e-n_b)+n_b.$$

$$2n-n_b+(2-n_{evf})n_f=2\chi ,$$

$$6(n-n_b)+3n_b=3n_f\Rightarrow n_f=2n-n_b$$

$$4(n-n_b)+2n_b=4n_f\Rightarrow n_f=n-\frac{n_b}{2}$$

$$4n_f=4(n-n_b-n_k)+2(n_b-n_{bk})+(4-k)n_k+(2-k)n_{bk},$$

$$\chi =\frac{k}{4}(n_k+n_{bk}),$$

$$\chi =\sum _k\frac{k}{4}(n_k+n_{bk})=\sum _{i=1}^N\text{index}({\mathbf{x}}_i).$$

$${\theta }_f(\mathbf{x})=\mathrm{d}{\theta }_f({\mathbf{x}}_a,\mathbf{x})$$

$$\underset{\mathbf{x}\to {\mathbf{x}}_b}{\lim }{\theta }_f(\mathbf{x})={\theta }_f({\mathbf{x}}_b)$$

$${\mathbf{e}}_1=\frac{{\mathbf{x}}_a-{\mathbf{x}}_i}{|{\mathbf{x}}_a-{\mathbf{x}}_i|},{\mathbf{e}}_2=\mathbf{n}\times {\mathbf{e}}_1,$$

$$r({\mathbf{x}})=|{\mathbf{x}}-{\mathbf{x}}_{i}|\quad,\quad\theta({\mathbf{x}})=\mathrm{atan2}\Big{(}({\mathbf{x}}-{\mathbf{x}}_{i})\cdot{\mathbf{e}}_{2},({\mathbf{x}}-{\mathbf{x}}_{i})\cdot{\mathbf{e}}_{1}\Big{)}$$

$$f:S^1\mapsto S^1/Q.$$

$$\mathrm{index}({\mathbf{x}}_i)=\frac{1}{2\pi }{\oint }_{{\mathcal{C}}_i}\mathrm{d}{\theta }_f$$

$$\mathrm{index}({\mathbf{x}}_i)=\frac{1}{2\pi }\sum _{p=1}^{v_i}\left({\alpha }_p-\frac{\pi }{2}\right)=\frac{1}{2\pi }(2\pi -v_i\frac{\pi }{2})=\frac{4-v_i}{4},$$

$$\sum _{i=1}^N\text{index}({\mathbf{x}}_i)=\chi$$

$$\mathbf{f}=(\cos 4{\theta }_f,\sin 4{\theta }_f)\equiv (f_1,f_2).$$

$$f=f_1+i\text{leavevmode}{f}_2$$

$$\mathbf{f}=(\cos 4{\theta }_n,\sin 4{\theta }_n)\mathrm{on}\Gamma .$$

$${\theta }_f=\frac{atan2(f_2,f_1)}{4}.$$

$$\sum _{i=1}^N\text{index}({\mathbf{x}}_i)=4\chi .$$

$$E(f_1,f_2)=\underset{\text{smoothing}}{\underbrace{\frac{1}{2}{\int }_{\mathcal{S}}\left({\left|\nabla f_1\right|}^2+{\left|\nabla f_2\right|}^2\right)d\mathcal{S}}}+\underset{\text{penality}}{\underbrace{\frac{1}{4{ϵ}^2}{\int }_{\mathcal{S}}{\left(f_1^2+f_2^2-1\right)}^{2}d\mathcal{S}}}.$$

$${\nabla }^2{f}_i-\frac{1}{{ϵ}^2}\left(f_1^2+f_2^2-1\right)f_i=0i=1,2.$$

$$E=\pi \left(\sum _{i=1}^N\text{index}({\mathbf{x}}_i^c{)}^2\right)\log (1/ϵ)+W+\mathcal{O}(1/|\log ϵ|).$$

$$W=-\pi \sum _{i=1}^N\sum _{\underset{j\ne i}{j=1}}^N\text{index}({\mathbf{x}}_i^c)\text{leavevmode}\text{index}({\mathbf{x}}_j^c)\log |{\mathbf{x}}_i^c-{\mathbf{x}}_j|+\mathcal{R}$$

$$\left\{{\mathbf{e}}^p={\mathbf{p}}_j-{\mathbf{p}}_i,{\mathbf{t}}^p={\mathbf{n}}^p\times {\mathbf{e}}^p\right\}$$

$$f_1^p=\cos 4{\theta }_f^p,f_2^p=\sin 4{\theta }_f^p,$$

$${\omega }^{(1)}(\xi ,\eta )=1-2\eta ,{\omega }^{(2)}(\xi ,\eta )=2(\xi +\eta )-1,{\omega }^{(3)}(\xi ,\eta )=1-2\xi ,$$

$${\Theta }^{(1)}={\theta }_f^{(1)},{\Theta }^{(2)}={\theta }_f^{(2)}+{\alpha }^{(2)},{\Theta }^{(3)}={\theta }_f^{(3)}+{\alpha }^{(3)}.$$

Full text

1 Introduction

The Finite element method (FEM) provides a powerful and versatile framework for numerical simulation, which however heavily relies on mesh generation, the decomposition of a geometrical region into simple shaped finite elements. In two-dimensional geometries, two kinds of elements exist: triangles and quadrangles. Quadrangular meshes are deemed better than triangular ones because ($i$) there are half as many quadrangles than triangles for the same number of vertices; ($ii$) it is possible to define tensorial operations on quadrangles; and ($iii$) quadrangular meshes ease the tracking of preferred directions in mesh refinement.

However, the generation of quadrangular meshes remains a challenging task, for which many strategies have been explored. Some of them, based on surface parameterization, are suitable for the generation of structured quadrangular meshes, close to regular (square) grids. A crossfield may be used to determine the appropriate parameterization, either on a patch [1] or globally [2]. A crossfield can also be used for partitioning the surface into a set of curvilinear quadrangular regions (a polyquad), then trivially quadrangulable [3]. The parameterization can also be deduced from a singularity graph [4]. In this paper, the primary concern is however to use crossfields as part of another meshing strategy: a frontal approach firstly proposed by [5] that consists in recombining triangles into quadrangles. This can be done efficiently [6] but the quality of the quadrangles strongly depends on the node location. A heuristic to obtain well distributed nodes is to spawn them following consistent directions, such as those suggested by a smooth crossfield. Such a frontal approach allows building unstructured quadrangular meshes with varying element size. Other advantages of quadrilateral meshes exist for specific finite element models: for examples, triangular plate bending elements are stiffer than quadrilateral ones with the same number of vertices

Although there exist various ways to represent discrete crossfields [7, §5], their computation generally relies on some smoothing process, possibly under constraints. For an angle-based representation, a crossfield is pictured as four orthogonal or opposite vectors. From this representation, it is possible to formulate the quadrangulation as a mixed-integer problem [1]. More advanced mathematical notions such as holonomy [8] may be used as well to design crossfields. This approach requires to build a metric on the 2-manifold.

In this paper, the so-called Cartesian (complex) representation [9] is adopted. This representation naturally takes the symmetries of the cross into account, and the crossfield is identified with a complex-valued function. Complex analysis gives then a large and useful background, especially about the theoretical analysis of critical points. The second term of the Ginzburg-Landau functional is controlled by a parameter depending on the local mesh size. When this parameter is small enough, the minimization of the functional results in a smooth crossfield whose critical points are optimally located and whose critical points have indices with minimal absolute values, according to the theory. The previous approach closest to ours is that in [3]. It has only the energy term, but the vector field is constrained to have a norm close to the unity. Critical points are identified in this approach by computing an argument (angle) from the vector field, whereas we only need to compute the vector field norm, critical points being in our approach points where the crossfield norm locally vanishes.

Our main contribution is to express the crossfield problem with Ginzburg-Landau equations. Those equations rely on an interesting mathematical and physical backgrounds. In order to grasp the great understanding that Ginzburg-Landau functional provides to the crossfield problem, we first recall the topological constraints of full quadrangular (and triangular) mesh in section 2 and the link with cross (and asterisk, respectively) field in section 3. In section 4, we develop the intuition of using the Ginzburg-Landau functional for the crossfield problem and we give the related Ginzburg-Landau theory. We derive in section 5 a simple FEM scheme from the Ginzburg-Landau equations. Our numerical scheme is validated on the unit disk in regards with Ginzburg-Landau theory, section 6. On the 2-sphere section 7 we get a surprising but correct result. In section 8, the Ginzburg-Landau equations are modified to get better results on NACA profiles. Finally, we apply our simple finite scheme on the torus in section 9.

2 Topology of Triangular and Quadrilateral Meshes

Assume an orientable surface $\mathcal{S}$ embedded in $\Re^{3}$. Let $g$ be the number of handles of the surface. The topological characteristic $g$, which is also called the genus of the surface, is the maximum number of cuttings along non-intersecting closed curves that won’t make the surface disconnected. Let also $b$ be the number of connected components of the boundary $\partial\mathcal{S}$ of the surface. The Euler characteristic of $\mathcal{S}$ is then the integer

[TABLE]

One has $\chi=2$ for a sphere, whereas $\chi=1$ for a disk ($b=1$), and $\chi=0$ for a torus ($g=1$) or a cylinder ($b=2$).

Consider now a mesh on $\mathcal{S}$ with $n$ nodes (also called vertices), $n_{e}$ edges and $n_{f}$ facets. The Euler formula

[TABLE]

provides a general relationship betweeen the numbers of nodes, edges and facets in the mesh (details in [10]). If $n_{b}$ nodes (and hence $n_{b}$ edges) are on the boundary $\partial\mathcal{S}$, and if the number of edges (or nodes) per facet is noted $n_{evf}$ ($n_{evf}=3$ for triangulations and $n_{evf}=4$ for quadrangulations, meshes mixing triangles with quadrangles being excluded), the following identity holds : all facets have $n_{evf}$ edges, $n_{e}-n_{b}$ edges have two adjacent facets and $n_{b}$ edges have one adjacent facet. Hence the relationship

[TABLE]

Elimination of $n_{e}$ between (2) and (1) yields

[TABLE]

which is true for any triangulation or quadrangulation.

A regular mesh has only regular vertices. An internal vertex is regular if it has exactly $6$ adjacent triangles or $4$ adjacent quadrangles, whereas a boundary vertex is regular if it has exactly $3$ adjacent triangles or $2$ adjacent quadrangles. One has then

[TABLE]

and

[TABLE]

respectively for a regular triangulation and a regular quadrangulation, from a topological point of view. Substitution of (4) and (5) into (3) shows that only surfaces with a zero Euler characteristic can be paved with a regular mesh. If $\chi\neq 0$, irregular vertices will necessarily be present in the mesh.

The number and the index of the irregular vertices is tightly linked to the Euler characteristic $\chi$, which is a topological invariant of the surface. We call valence of a vertex the number of facets adjacent to the vertex in the mesh. In a regular mesh, all vertices have the same valence $v_{reg}$. In a non regular mesh, on the other hand, a number of irregular vertices have a valence $v\neq v_{reg}$, and one notes the integer $k=v_{reg}-v$ the valence mismatch of a vertex.

Assume a quadrangulation with $n_{k}$ irregular internal vertices of valence $v=4-k$, and $n_{bk}$ irregular boundary vertices of valence $2-k$, $k$ given. All other vertices are regular. There are then $n-n_{b}-n_{k}$ regular internal vertices of valence $4$, and $n_{b}-n_{bk}$ regular boundary vertices of valence $2$, so that one can write

[TABLE]

and the substraction of (3) with $n_{evf}=4$ yields

[TABLE]

showing that, in a quadrangulation, each irregular vertex counts for $\text{index}({\mathbf{x}}_{i})=k/4$ in the Euler characteristic, a quantity called the indice of the irregular vertex ${\mathbf{x}}_{i}$.

Summing up now on different possible values for $k$, one can establish that a quadrangulation of a surface with Euler characteristic $\chi$ verifies

[TABLE]

Consider, for instance, the quadrangulation of a disk, which is a surface with $\chi=1$. A minimum of $n_{1}=4$ irregular vertices of index $1/4$ must be present. They can be located either on the boundary (vertices of valence 1) or inside the disk (vertices of valence 3), Fig. 1.

Fig. 2 shows three different quadangulations of a L-shaped domain ($\chi=1$). Regular boundary nodes should all have a valence of 2. The mesh on the left has $6$ irregular vertices located at the corners of the domain : five with index $1/4$, and one with index $-1/4$. The central mesh, on the other hand, has the minimum amount of irregular vertices, i.e. four ones of index $1/4$. The right mesh generated by recombination of a standard Delaunay triangular mesh (see [6]) has twelve vertices of index $1/4$, and eight vertices of index $-1/4$, both on the boundary and inside the domain. Quality meshes should have as few irregular vertices as possible. In what follows, a general approach allowing to compute the position of such irregular vertices before meshing the surface is presented.

3 Why Crossfields?

Crossfields are auxiliary in the generation of quadrangular meshes. We shall show that nonregular vertices defined in the previous section are precisely the critical points of a crossfield, and that these critical points of the crossfield can also be related to the Euler characteristic of the meshed surface. This result represents an important theoretical limit on the regularity of quadrangular meshes.

3.1 Continuity

A crossfield $f$ is a field defined on a surface $\mathcal{S}$ with values in the quotient space $S^{1}/Q$, where $S^{1}$ is the circle group and $Q$ is the group of quadrilateral symmetry. Pictorially, it associates to each point of the surface $\mathcal{S}$, which has to be meshed, a cross made of four unit vectors that are orthogonal with each others in the tangent plane $T\mathcal{S}$ of the surface.

A surface $\mathcal{S}$ can be identified with its tangent space in any neighborhood $\sigma\subset\mathcal{S}$ that is sufficiently small to have curvature effects negligible. This local identification of the surface with a vector space endows it with a natural parallel transport rule, so that the angular differential $\mathrm{d}\theta_{f}({\mathbf{x}}_{a},{\mathbf{x}}_{b})$ can be defined as the minimal angle, with its sign, between the branches of $f({\mathbf{x}}_{a})$ and any of the branches of $f({\mathbf{x}}_{b})$ for any pair of points ${\mathbf{x}}_{a},{\mathbf{x}}_{b}\in\sigma$ where $f$ is defined, Fig. 3. Taking now as reference the cross $f({\mathbf{x}}_{a})$, an angular coordinate

[TABLE]

can be defined for crosses in $\sigma$. The crossfield $f$ is deemed continuous (regular) at ${\mathbf{x}}_{b}$ if the limit

[TABLE]

exists (i.e. is unique). It is then equal to $\theta_{f}({\mathbf{x}}_{b})$. Isolated points ${\mathbf{x}}_{i},\leavevmode\nobreak\ i=1\dots N$, of $\mathcal{S}$ where the limit (9) does not exist are called critical points or zeros of the crossfield.

3.2 Index and degree

Although defined locally, the notion of continuity gives unexpectedly valuable information about the topology of $\mathcal{S}$, which is a nonlocal concept. To see this, consider a crossfield $f$ defined on a quadrangular element delimited by four (possibly curvilinear) edges. Assume the crossfield is parallel to the four edges (i.e. one of the four branches of the cross is parallel to the tangent vector of the edge at each point of the edge, except the extremities) and prolongates smoothly inside the quadrangle. This field is discontinuous at corners where edges do not meet at right angle, but it is continuous everywhere else. Making the same construction for all elements of a quadrangular mesh, one obtains a crossfield $f$ topologically identified with the quadrangular mesh, and that is continuous everywhere except at the vertices of the mesh. This field has thus got isolated critical points at mesh vertices, but not all critical points have the same significance. Some critical points have a specific topological value, associated with the notion of index.

To introduce the notion of index, an angular coordinate needs to be defined for points in a neighborhood $\sigma_{i}$ of a critical point ${\mathbf{x}}_{i}$. Picking up an arbitrary regular point ${\mathbf{x}}_{a}\in\sigma_{i}$, ${\mathbf{x}}_{a}\neq{\mathbf{x}}_{i}$, the local unit vector basis

[TABLE]

with ${\mathbf{n}}$ the normal to $\mathcal{S}$, is constructed, and hence a local polar coordinate system

[TABLE]

can be defined for points in $\sigma_{i}$.

A circular curve $\mathcal{C}_{i}$ of infinitesimal radius centered around the vertex ${\mathbf{x}}_{i}$ is now considered. As the angles $\theta({\mathbf{x}})$ (10) and $\theta_{f}({\mathbf{x}})$ (8) are precisely the elements of the groups $S^{1}$ and $S^{1}/Q$, respectively, the crossfield on $\mathcal{C}_{i}$ can be regarded as a mapping

[TABLE]

The mapping is continuous, since $\mathcal{C}_{i}$ circles around the critical point ${\mathbf{x}}_{i}$, but it does not cross it. The index of $f$ at ${\mathbf{x}}_{i}$ is the degree of the mapping (11), i.e. the number of times the codomain wraps around the domain under the mapping. Its algebraic expression is easily expressed in terms of the angles $\theta$ and $\theta_{f}$ as

[TABLE]

where $2\pi$ is $\oint_{\mathcal{C}_{i}}\mathrm{d}\theta$. In case of a vertex ${\mathbf{x}}_{i}$ of valence $v_{i}$, i.e. a vertex adjacent to $v_{i}$ quadrangular elements, the integral evaluates as

[TABLE]

where the $\alpha_{p}$’s are the angles of the $v_{i}$ quadrangular elements adjacent to the considered vertex ${\mathbf{x}}_{i}$, and where the obvious relationship $\sum_{p=1}^{v_{i}}\alpha_{p}=2\pi$ has been used. The crossfield $f$ has index 0 at vertices adjacent to four quadrangular elements, whereas it has index $1/4$ ($-1/4$) at vertices adjacent to 3 (5, respectively) quadrangular elements meet, Fig. 4. As one sees, the index is a topological characteristic of the crossfield $f$ at the critical point ${\mathbf{x}}_{i}$. It does not depend on the choice of the curve $\mathcal{C}_{i}$, nor on the choice of an angular reference for the angles $\theta({\mathbf{x}})$ and $\theta_{f}({\mathbf{x}})$.

3.3 Poincaré-Hopf theorem

Equation (12) relates the index of the crossfield at a critical point ${\mathbf{x}}_{i}$ with one fourth of valence $k_{i}=4-v_{i}$ of the corresponding mesh vertex. This result can be combined with the algebraic topology result of previous section (7) that each internal irregular vertex of valence $k_{i}$ counts for $k_{i}/4$ in the Euler characteristic of the underlying surface. This yields the relationship

[TABLE]

for the critical points of a crossfield $f$ defined on a surface $\mathcal{S}$.

This is a generalization Poincaré-Hopf theorem, which states that the sum of the indices of the critical points of a vector field ${\mathbf{v}}$ defined on a surface $\mathcal{S}$ without boundary is equal to the Euler characteristic of the surface. This famous theorem draws an unexpected and profound link between two apparently distinct areas of mathematics, topology and analysis. Whereas vector fields have integer indices at critical points, crossfields have indices that are multiples of 1/4. Still the topological relationship (13) of Poincaré-Hopf holds in both cases. Actually, our developments reach same inferences as [11].

4 Crossfield Computation: the Planar Case

We introduce the representation of a crossfield by means of a vector field. From this representation, we derive the problem to solve that corresponds to minimize Ginzburg-Landau functional. Its asymptotic behavior provides suitable critical points, if any.

4.1 Vector representation of crossfields

Only scalar quantities can be compared at different points of a manifold. For the comparison or, more generally, for differential calculus with nonscalar quantities like crossfields, a parallel transport rule needs to be defined on the manifold. On a surface (two-manifold), this rule can take the form of a regular vector field which gives at each point the direction of the reference angle 0. Poincaré-Hopf theorem says that such a field does not exist in general, and in particular on manifolds whose Euler characteristic is not zero. The situation is however easier in the planar case. A global Cartesian coordinate frame can always be defined over the plane, and be used to evaluate the orientation of the crossfield. We shall therefore expose the crossfield computation method in the planar case, and then generalize to nonplanar surfaces, where we will have to deal with local reference frames, in a subsequent section.

A cross $f({\mathbf{x}})$ is an element of the group $S^{1}/Q$, which can be represented by the angle $\theta_{f}({\mathbf{x}})$ it forms with the local reference frame. Yet, due to the quadrilateral symmetry, four different angles in $[0,2\pi[$ represent the same crossfield $f({\mathbf{x}})$. Let for instance the angles $\theta_{1}=0$ and $\theta_{2}=\pi/2$ represent the same cross. The average $(\theta_{1}+\theta_{2})/2=\pi/4$ represents another cross, whereas the difference $\theta_{2}-\theta_{1}=\pi/2$ is not zero. So, we have $\frac{1}{2}(x+x)\neq x$ and $x-x\neq 0$, which clearly indicates that the values of the crossfield $f$ do not live in a linear (affine) space. This makes the representation by $\theta_{f}$ improper for finite element interpolation. The solution is two-fold. First, the angle $\theta_{f}$ is multiplied by four, so that the group $S^{1}/Q$ is mapped on the unit circle $S^{1}$, and the cross $f$ is therefore represented by a unit norm vector ${\mathbf{f}}$. Then, the vector is represented in components in the reference frame as

[TABLE]

This vector may be represented by a complex-valued function

[TABLE]

This representation corresponds to a vector field that is described by a complex exponential whose argument is $4\leavevmode\nobreak\ \theta$. A crossfield is thus depicted by the fourth roots of a (unit) complex number. This observation may be generalized for directional fields with $n$ symmetries [7, §5.2].

4.2 Laplacian smoothing

Computing the crossfield $f$ consists thus now of computing the vector field representation ${\mathbf{f}}$, which obviously lives in a linear space (a 2D plane). The components of ${\mathbf{f}}$ are fixed on the boundaries of $\Gamma=\partial\mathcal{S}$ so that the crosses are parallel with the exterior normal vector ${\mathbf{n}}=(\cos\theta_{n},\sin\theta_{n})$ i.e.

[TABLE]

Propagating ${\mathbf{f}}$ inside $\mathcal{S}$ is here done by solving a Laplacian problem. Even though the vector representation ${\mathbf{f}}$ is unitary on $\Gamma$, it tends to drift away from $S^{1}$ inside the domain. The computed finite element solution ${\mathbf{f}}$ lies therefore outside the unit circle and must be projected back on $S^{1}$ to recover the angle

[TABLE]

Due to the multiplication by 4, the indices of the critical points of the vector field ${\mathbf{f}}$ verify

[TABLE]

4.3 The Ginzburg-Landau model

Numerical experiments show that the norm of the vector field ${\mathbf{f}}$ computed by Laplacian smoothing (see previous section) decreases quite rapidly as one moves away from the boundary $\partial S$, leaving in practice large zones in the bulk of the computational domain where the solution is small, and the computed crossfield inaccurate, Fig. 5a. A more satisfactory formulation consists of ensuring that the norm of ${\mathbf{f}}$ remains unitary over the whole computational domain, Fig. 5b. This problem can be formulated in variational form in terms of the Ginzburg-Landau functional

[TABLE]

The first term minimizes the gradient of the crossfield and is therefore responsible for the laplacian smoothing introduced in the previous section. The second term is a penality term that vanishes when ${\mathbf{f}}\in S^{1}$. The penality parameter $\epsilon$, called coherence length, has the dimension of a length. The Euler-Lagrange equations of the functional (15) are the quasi-linear PDE’s

[TABLE]

called Ginzburg-Landau equations. If $\epsilon$ is small (enough) with respect to the dimension of $\mathcal{S}$, then ${\mathbf{f}}$ is of norm $1$ everywhere but in the vicinity of the isolated critical points ${\mathbf{x}}_{i}^{c}$.

The asymptotic behavior of Ginzburg-Landau energy can be written as

[TABLE]

with

[TABLE]

as $\epsilon\rightarrow 0$ (see [12], Introduction, Formulae 11 and 12).

In asymptotic regime, the energy is thus composed of three terms. The first term of (17) blows up as $\epsilon\rightarrow 0$, i.e. energy becomes unbounded if critical points are present. When $\epsilon$ is small, this first term dominates, and one is essentially minimizing $\sum_{i=1}^{N}\text{index}({\mathbf{x}}_{i}^{c})^{2}$ with the constraint (14). This indicates that a critical point of index $2$ has a cost of $4$ in terms of energy, whereas 2 critical points of index $1$ have a cost of $2$. All critical points should therefore be of index $\pm 1$, and their number should be $N=4\leavevmode\nobreak\ |\chi|$. This is indeed good news for our purpose : good crossfields should have few critical points of lower indices.

The second term of (17) is the renormalized energy $W$ (18). It remains bounded when $\epsilon$ tends to [math]. The double sum in $W$ reveals the existence of a logarithmic force between critical points. The force is attractive between critical points with indices of opposite signs, and repulsive between critical points with indices of the same signs. The second term in (18) is more complicated and is detailed in [12]. Basically, $\mathcal{R}$ represents a repulsing force that forbids critical points to approach the boundaries.

Finally, the third term in (17) vanishes as $\epsilon\rightarrow 0$. At the limit, all energy is thus carried by the critical points of the field. All this together allows to believe that Ginzburg-Landau model is a good choice for computing crossfields. It produces few critical points, which are moreover well-distributed over the domain.

5 Computation of Crossfields: Nonplanar Generalization

The finite element computation method for crossfields is now generalized to the case of nonplanar surfaces. Consider the conformal triangulation $\mathcal{S}=\cup_{ijk}\Omega_{ijk}$ of a nonplanar surface manifold $\mathcal{S}$, each triangle $\Omega_{ijk}$ being defined by the vertices ${\mathbf{p}}_{i}$, ${\mathbf{p}}_{j}$ and ${\mathbf{p}}_{k}$. Since no global reference frame exists on a nonplanar surface, a local reference frame is associated to each edge of the triangulation. Let $e^{p}$ be the $p^{\mathrm{th}}$ edge of the mesh, joining nodes ${\mathbf{p}}_{i}$ and ${\mathbf{p}}_{j}$, and ${\mathbf{n}}^{p}$ be the average of the normals vectors of the two triangles adjacent to $e^{p}$. The vectors

[TABLE]

form a local frame $\{\hat{e}^{p},\hat{t}^{p}\}$ which enables the representation of the connector values of the discretized crossfield ${\mathbf{f}}$,

[TABLE]

which are attached to the center of the edges of the triangulation. Actually, $\theta_{f}^{p}$ is assumed to be the same along $\bm{e}^{p}$ within both planes of triangles sharing $e^{p}$. This assumption eases computations and gives a planar-like representation, Fig. 6a.

As the connector values are attached to the edges of the mesh, and not to the nodes, Crouzeix-Raviart interpolation functions are used instead of conventional Lagrange shape functions, [13]. The Crouzeix-Raviart shape functions $\omega^{p}$ equal $1$ on corresponding edge $e^{p}$, and $-1$ on the opposite vertices (Fig. 7) in the two adjacent triangular elements. They are polynomial and their analytic expression in the reference triangle $\{\xi\in[0,1],\leavevmode\nobreak\ \leavevmode\nobreak\ \eta\in[0,1-\xi]\}$ reads

[TABLE]

where indices $(1),(2)$ and $(3)$ enclosed in parentheses denote the local edge numbering in the considered triangular element.

Each of the three edges of a triangle $\Omega_{ijk}$ has its own local reference frame. If one is to interpolate expressions involving the vector field ${\mathbf{f}}$ over this element, the three edge-based reference frames have to be appropriately related with each other (see [14]). We arbitrarily take the reference frame of the first edge of the element as reference, and express the angular coordinate of the two other edges in function of this one with the relationships (Fig. 6b)

[TABLE]

Thus, the 6 local unknowns of triangle $\Omega_{ijk}$ can be expressed as a function of the 6 edge unknowns by

[TABLE]

and we have the interpolation

[TABLE]

for the vector field ${\mathbf{f}}$ in the triangle $\Omega_{ijk}$.

A Newton scheme is proposed to converge to the solution. The Newton iteration at stage $n$ for solving (16) consists of solving:

[TABLE]

The $6\times 6$ elementary matrix ${\mathbf{K}}^{ijk}$ and the $6\times 1$ elementary vector ${\mathbf{B}}^{ijk}$ of element $\Omega_{ijk}$ are then given by ${\mathbf{K}}^{ijk}$

[TABLE]

and

[TABLE]

with $m,n=1\dots 3$.

It is then necessary to transform those elementary matrix and vector in the reference frames of the edges as

[TABLE]

Then, standard finite element assembly can be performed. Boundary conditions are simply

[TABLE]

on every edge $e^{p}$ of $\partial\mathcal{S}$. This nice simplification is due to the fact that unknowns are defined on the reference frame of the edges.

6 Numerical Validation: the Unit Disk

We compute the analytical location of critical points of a directional field defined on the unit disk. The calculations are based on the Ginzburg-Landau results, described in section 4.3. The numerical location obtained by our FEM is compared to the analytical one.

Let $\mathcal{S}$ be the open unit disk in $\Re^{2}$, i.e.

[TABLE]

For a star-shaped planar domain such as $\mathcal{S}$ with a smooth boundary $\partial\mathcal{S}$ of exterior normal $\nu$ and tangent $\tau$, whose vector field has $d$ critical points of index $+1$ at $X^{c}=\{\mathbf{x}_{1}^{c},...,\mathbf{x}_{d}^{c}\}\in\mathcal{S}$, the asymptotic energy $E_{\epsilon}$ (in complex form) becomes

[TABLE]

where $W(X^{c})$ is the renormalized energy

[TABLE]

where $\Phi$ is given by the following Neumann problem

[TABLE]

and $R$ is the regular part of $\Phi$:

[TABLE]

$E_{\epsilon}$ is minimum when the critical points are located appropriately, i.e. when (23) is minimum. The renormalized energy $W$ corresponds to the Ginzburg-Landau energy (22) when the singular core energy $\pi d|\log(\epsilon)|$ has been removed. Since $W$ depends only on the location of the critical points, it is possible to compute their location in the case of the unit disk, in order to get an optimal directional field.

The minimum of $W$ is obtained by sampling points within the unit disk. It is assumed that the $d$ critical points exhibit the $d$ symmetries of their group (the quadrilateral group in the case $d=4$). In other words, it means that they are at the same distance $r^{c}$ from the center of the disk (i.e. the origin $(0;0)$), and separated two-by-two with an angle of $2\pi/d$ radians.

The Neumann problem (24) is solved by decomposing $\Phi=\Phi^{0}+\Phi^{1}$. The first term $\Phi^{0}$ is the Green function of a two-dimensional Laplacian operator, while the second one $\Phi^{1}$ is obtained by separation of variables $(r,\theta)$. The solution is then

[TABLE]

where $A_{i,n}$ depends on the location of the i-th critical point, which is parameterized by $r^{c}$. It is possible to show that the second term of (23) is zero, Appendix B. The analytical solution of Neumann problem is derived into the Appendix A.

The evaluation of $W$ consists of computing the first and last terms, by sampling the disk. The sampling is done by selecting $d$ critical points spaced by $2\pi/d$ radians. The distance $r^{c}$ is sampled between zero and one. The distance $r^{c^{*}}$ which gives the lowest value of $W$ defines the location of the critical points. A Python script performs the evaluations and returns the optimal distance $r^{c^{*}}$, Fig. 8.

The corresponding directional fields are computed, and their critical point locations are compared with circles which radii correspond to $r^{c^{*}}$, Fig. 9. The location of critical points are really close to the estimation based on the analytical solution of $W(X^{c})$ in the case of the unit circle. They tend to draw the corners of the polygon of symmetry: a square in the case of the crossfield, Fig. 9a and a regular hexagon for the asterisk field, Fig. 9b. The critical points are quite close to the unit circle. The more critical points, the closer to the unit circle they are. We understand that the repulsion term is stronger than the regularization term within the domain. The regularization term is only able to forbid critical points to be on the boundary, i.e. the unit circle.

7 A Surprising Result: the Sphere

Let us compute the crossfield on a unit sphere. The sphere has no boundary so we choose randomly one edge of the mesh and fix the crossfield for this specific edge. The mesh of the sphere is made of 2960 triangles (see Fig. 10). A value of $\epsilon=0.1$ was chosen for the computation. A total of $29$ Newton iterations were necessary to converge, by reducing the residual norm to $10^{-12}$.

The location of the $8$ critical points is indeed not what we expected: our initial intuition was that critical points would be located at the corners of an inscribed cube of side $1/\sqrt{3}$. In all our computations i.e. while changing the mesh and $\epsilon$, critical points are located on two squares of side $1/\sqrt{3}$, those two squares being tilded by $45$ degrees around their common axe (see Fig. 10). Equilateral triangle patterns are formed between critical points that belong to both squares. In reality, our solution is the right solution. In the asymptotic regime, the location ${\mathbf{x}}_{i}^{c}$ of the $8$ critical points tends to minimize $-\sum_{i}\sum_{i\neq j}\log|{\mathbf{x}}_{i}^{c}-{\mathbf{x}}_{j}^{c}|$ (see Equations (17) and (18)). We have thus computed $-\sum_{i}\sum_{i\neq j}\log|{\mathbf{x}}_{i}^{c}-{\mathbf{x}}_{j}^{c}|$ for tilting angles ranging from [math] to $\pi/2$. Fig. 11 shows clearly that the minimum of the energy corresponds to an angle of $\pi/4$, which is exactly what is found by the finite element formulation.

Fig. 12 shows the crossfield as well as the separatrices. The separatrices were computed “by hand”.

The solution that has been found is related to what is called the Whyte’s problem (cf. [15, 16]) that consists of finding $N$ points on the sphere which positions maximize the product of their distances. The critical points are called logarithmic extreme points or elliptic Fekete points (see [17]).

The specific configuration that corresponds to $N=8$ is called an anticube (or square antiprism) and is exactly the one that was found numerically.

In the case of an asterisk field, the critical points are the summits of an icosahedron, which is the solution of Whyte’s problem for $N=12$. This superb result shows that it is indeed possible to use crossfields not only for building quadrangles but also to build equilateral triangles.

Actually, it is possible to show that the critical points computed over the sphere by Ginzburg-Landau correspond to the solution of Whyte’s problem for any even value of $N$ (see [18]).

8 Weak Boundary Conditions

In this section, we have computed the graph of singularities of a standard CFD test case: a three component wing domain with $\chi=-2$. This example is very similar to the one presented by [3]. The solution has been computed on a non uniform triangular mesh of about $15,000$ triangles. The graph of singularities has been depicted on Fig. 14.

Weak boundary conditions have been applied to the different components of the wing where a penalization replaces the strong imposition of $f$ on boundaries. A new term is thus added to Energy (15) for taking into account boundary conditions:

[TABLE]

where $\bar{f}_{1}$ and $\bar{f}_{2}$ are values of the crosses that are weakly imposed on the boundary and $L$ the characteristic size of the problem. This new treatment allows singularities to migrate on the boundary, making their repulsive action finite. Figure 14 clearly shows that effect: a singularity of index $1/4$ sits on the leading edge of the slat, allowing a clean decomposition of the domain. The same migration is also observed on the leading edge of the profile. A strong imposition of boundary conditions naturally leads to singularities that are very close to regions of the boundary with high curvature, usually at a distance from the boundary that is one mesh size. Artificial boundary layers are thus added to the decomposition (see [3, Fig. 12 and 14]).

9 Application of our FEM Scheme to the Torus

The Euler-Poincaré characteristic of the torus is $\chi=0$. Theoretically, we should obtain a crossfield without critical points. But our FEM scheme gives crossfield with twelve critical points, located where the Gaussian curvature is maximal (exterior) or minimal (interior), Fig. 15. Fig. 15a shows that the six critical points located on the maximal Gaussian curvature line are facing the six corresponding critical points located on the minimal Gaussian curvature line. Moreover, as the former have an index $+1/4$, and the latter an index $-1/4$, Fig. 15b, the index sum of the surface is zero, as predicted by the Poincaré-Hopf theorem.

Our FEM scheme does not reach however the asymptotic behavior $(\epsilon\rightarrow 0)$ of the Ginzburg-Landau functional. It means that our penalty factor $\epsilon$ is not low enough. Otherwise, the computed crossfield should not have any critical points owing to (17). Actually, the computed crossfield has a lower energy ($72.10$) than the crossfield with no critical point that could be drawn by aligning crosses with the main curvatures of the surface ($84.58$). The tentative polyquad decomposition shown in Fig. 15c indicates that the field computed with the Ginzburg-Landau approach tends to be more uniform, in order to reduce the Dirichlet energy. It confirms that the Dirichlet term is stronger than the penalty term.

10 Conclusion

This article has demonstrated the consistency of the Ginzburg-Landau theory to compute directional fields on arbitrary surfaces. The proposed approach relies on a physical and mathematical backgrounds. This provides proofs, analytical solutions and helps delineating fundamental mathematical properties that can be exploited in algorithms.

In particular, the Ginzburg-Landau theory states that when the coherence length $\epsilon$ is small enough, the asymptotic behavior is reached, i.e., the number of critical points of the crossfield is minimal, their index is also minimal and they are optimally distributed. A simple FEM scheme has been implemented to validate numerically this assertion. Crossfields have been computed on the unit disk and solutions conform with the Ginzburg-Landau theory have been found. The location of computed critical points on the 2-sphere corresponds to the solution of Whyte’s problem: for a crossfield they are at the summits of an anticube whereas for an asterisk field they are at the summits of a regular dodecahedron.

By weakening the boundary conditions of the Ginzburg-Landau problem, critical points are no longer repelled in the interior of the domain and can be located on the boundary, which improves the polyquad decomposition in the case of the NACA profiles.

Finally, the process is applied to the quadrangular meshing of the coastal domain around Florida peninsula, Fig. 16. Quadrangles are merged from right-angled triangles whose vertices have been spawned along the integral lines of a crossfield, Fig. 17a. One sees on Fig. 17b how the edges of the recombined quadrangular elements tend to follow the crossfield, and the final mesh is of satisfying quality, Fig. 18.

The input triangular mesh can be improved by using an asterisk field. This field is used to spawn vertices which are consistent with an equilateral triangular grid, Fig. 19a. The vertices tend to have the correct valence, except in some regions where the size field changes, Fig. 19b. The final triangular mesh exhibits a smoother distribution of equilateral triangles through the domain, while the mean quality $\bar{\gamma}$ has been improved to 0.9559 (from $\bar{\gamma}=$0.9505 for the initial mesh), Fig. 20.

Further work will focus on highly improving the numerical scheme that solves Ginzburg-Landau equations, in order to make it competitive.

Acknowledgements

The present study was carried out in the framework of the project ”Large Scale Simulation of Waves in Complex Media”, which is funded by the Communauté Française de Belgique under contract ARC WAVES 15/19-03 with the aim of developing and using Gmsh ([19]).

Appendix

This appendix provides an analytical form of the renormalized energy $W(X^{c})$ of equation (24) for a unit disk . We first compute $\Phi$ and then prove that second term of (23) is equal to zero in the case of a unit circle.

Appendix A Solving the Neumann Problem i.e. Computing $\Phi(x)$ of Equation (24)

Assume a unit circle $S$. The analytical value of $f$ on the boundary $\partial S$ of $S$ is

[TABLE]

as one direction has to be aligned with $\tau$ along the circle. The Neumann boundary condition is thus

[TABLE]

since $f\times\nabla f\cdot\tau=d$ on $\partial\mathcal{S}$. Indeed, from

[TABLE]

and

[TABLE]

the condition (28) corresponds to the imaginary part of the corresponding complex product.

The Green function of the two-dimensional Laplacian operator is part of $\Phi$

[TABLE]

Even if $\nabla^{2}\Phi^{0}=2\pi\sum_{i=1}^{d}\delta(\mathbf{x}-\mathbf{x}_{i}^{c})$, the flux (per unit of length) $\nabla\Phi^{0}\cdot\nu$ does not correspond to (28). The solution $\Phi$ contains another term $\Phi^{1}$. It may also be written as a sum of the contributions coming from the $d$ critical points. Therefore,

[TABLE]

such that

[TABLE]

Function $\Phi^{1}_{i}$ can be written as series of circular harmonics

[TABLE]

where $(r;\theta)$ are polar coordinates. We search for the solution of a Neumann problem which is defined to a constant. Setting set $A_{i,0}=0$ assigns to zero the average of $\Phi^{1}_{i}$. The idea is simple. We use $\Phi^{1}$ which is harmonic to remove all oscillatory parts of $\nabla\Phi\cdot\nu$ along the boundary $\partial\mathcal{S}$.

Let us assume the i-th critical point is located on the x axis (the real axis), i.e. $\mathbf{x}_{i}^{c}=(r^{c};0)$ with cartesian coordinates (see Figure 21).

One has

[TABLE]

and

[TABLE]

The last expression can be reformulated as

[TABLE]

with $s=2r^{c}/(1+{r^{c}}^{2})$. Taking into account the identity

[TABLE]

we have

[TABLE]

and

[TABLE]

Powers of $\cos(\theta)$ appear in (31). In order to replace such powers by $\cos(n\leavevmode\nobreak\ \theta)$’s like in Equation (32), we use a well known property of Chebyshev polynomials : $P_{n}(\cos(\theta))=\cos(n\leavevmode\nobreak\ \theta)$. We thus have

[TABLE]

where the ${\mathcal{P}}_{mn}$’s are the entries of the Chebyshev coefficient matrix $\mathcal{P}$. Equation (33) can thus be regarded as a system of equations

[TABLE]

The system matrix $\mathcal{P}$ is lower triangular, so the system can be inverted easily

[TABLE]

or equivalently, back with the initial notation,

[TABLE]

Thus,

[TABLE]

with

[TABLE]

Finally, we get the following series for the normal derivative of $\Phi_{i}^{0}$:

[TABLE]

We get the final condition

[TABLE]

The boundary condition should be non oscillatory: So,

[TABLE]

Finally

[TABLE]

Appendix B $\displaystyle{\int_{\partial\mathcal{S}}\Phi\leavevmode\nobreak\ f\times f_{\tau}\leavevmode\nobreak\ ds}$ is Zero along a Circle

We want to show that

[TABLE]

when $\partial\mathcal{S}$ is a unit circle.

In the previous section, we have shown that

[TABLE]

Besides, $\Phi^{1}$ has been derived such that it is non oscillatory along the unit circle. Hence, it remains to show

[TABLE]

We can express that integral with complex variables

[TABLE]

with the complex logarithm, which has two features:

• •

the complex logarithm is a multivalued function

• •

the complex logarithm has a peculiar singularity in zero

Those features are due to the fact that zero is a branch point. In our case, the branch points are the critical points $z_{k}^{c}$. A branch cut has to be drawn for each critical point. If $z_{k}^{c}=r^{c}\exp(i\leavevmode\nobreak\ \gamma_{k})$, the branch cut is such that $\theta\in[\gamma_{k};\gamma_{k}+2\pi),z=r\exp(i\leavevmode\nobreak\ \theta)\in\mathbb{C}$ (red line on Fig. 22). The logarithm is analytical thanks to the branch cut, and its contour integral may be written as

[TABLE]

The antiderivative of $\log(w)$ is $w\left(\log(w)-1\right)$, yielding the following

[TABLE]

The last step is due to the fact that $\gamma_{k}^{-}-\gamma_{k}^{+}=2\pi$ (the branch cut, in other words).

[TABLE]

which depends on the branch cut.

The $d$ critical points $z_{k}^{c}$ being at the same distance $r^{c}$ from the origin and evenly spaced by an angles $2\pi/d$, we get

[TABLE]

which is zero since the sum of $d>1$ complex numbers corresponding to points evenly distributed along on a circle centered at the origin is zero.