Ginzburg-Landau energy and placement of singularities in generated cross
fields

Alexis Macq1,2, Maxence Reberol1, François Henrotte1, Pierre-Alexandre
Beaufort1,3, Alexandre Chemin1, Jean-François Remacle1, Jean Van Schaftingen1

1 UCLouvain, Louvain-la-Neuve, Belgium

2 Haute Ecole Galilée - ISPG, Brussels, Belgium

3 University of Bern, Switzerland

Contact<EMAIL_ADDRESS>
## Abstract

Cross field generation is often used as the basis for the construction of
block-structured quadrangular meshes, and the field singularities have a key
impact on the structure of the resulting meshes. In this paper, we extend
Ginzburg-Landau cross field generation methods with a new formulation that
allows a user to impose inner singularities. The cross field is computed via
the optimization of a linear objective function with localized quadratic
constraints. This method consists in fixing singularities in small holes
drilled in the computational domain with specific degree conditions on their
boundaries, which leads to non-singular cross fields on the drilled domain. We
also propose a way to calculate the Ginzburg-Landau energy of these cross
fields on the perforated domain by solving a Neumann linear problem. This
energy converges to the energy of the Ginzburg-Landau functional as epsilon
and the radius of the holes tend to zero. To obtain insights concerning the
sum of the inner singularity degrees, we give: (i) an extension of the
Ginzburg-Landau energy to the piecewise smooth domain allowing to identify the
positions and degrees of the boundary singularities, and (ii) an
interpretation of the Poincaré-Hopf theorem focusing on internal
singularities.

## Introduction

Automatic block-structured quadrilateral meshing techniques have made
significant progress in the last decade [1, 2, 3] thanks to broad use of
smooth cross fields [4, 5, 6] aligned with the model boundaries. The structure
of the resulting quadrilateral meshes is mainly determined by the cross field
properties, for instance, the mesh irregular vertices match the cross field
singularities and the quad sizes are related to the size map inherent to the
cross field. In this context, understanding and controlling the cross field
properties is the natural way to influence the quad mesh and to tailor it for
engineering applications. The optimization of the mesh characteristics (size
map, number of structured blocks, position of irregular vertices) via the
cross field is still largely unexplored but of crucial importance for
numerical simulation.

Throughout this paper, a cross field will be considered as a discrete object
that assigns four orthogonal directions and a norm to each point of a
discretized domain. Those fields may include singularities which are points
where four orthogonal directions can not be defined while keeping the field
smooth. Our paper focuses on the management of the singularities which is
primordial to adapt to the needs of a high-precision numerical simulation.
Indeed, the singularities appear to coincide with the corners of the blocks of
the subsequent block-structured quadrangular mesh (for such results linked
with Ginzburg-Landau theory see [7, 8]). The aim of this paper is to let a
user impose inner singularities in a cross field while giving him methods to
evaluate his choices in terms of Ginzburg-Landau energy.

This paper begins with a short literature review 1. We then propose
mathematical background 2 in which we give a problem to evaluate the choice of
singularities as perforations in the computational domain 2.6. The asymptotic
equivalence between this energy and the Ginzburg-Landau energy functional is
given. We also highlight a strong link between the energetic Ginzburg-Landau
functional and an energy that involves only some positions and degrees.
Indeed, the positions and degrees that minimize this last energy are the same
as the positions and degrees of the singularities of the minimizers of the
Ginzburg-Landau functional 2.6.1. Next, we introduce the main proposed method
of this article 3.1. It consists in trapping the singularities in small holes
being drilled in the computational domain and holding specific boundary
conditions. By doing so, a user can choose a configuration of inner
singularities as holes in a computed non-singular smooth cross field.
Unfortunately, our main method does not allow to impose boundary
singularities. This have lead us to propose some theoretical novelties. They
bring clues to understand which imposed inner singularity configurations will
lead to respect the Poincaré-Hopf theorem in case of angular boundaries. In
terms of theoretical inputs, we extend the Ginzburg-Landau energy to
piecewise-smooth closed domains and give a strategy for managing singularity
configurations based on the minimization of this new Ginzburg-Landau boundary
energy 3.2. We also give new interpretation of the generalized Poincaré-Hopf
theorem focusing on internal singularities 4.2. As further thoughts, we
present a Neumann problem with zero-radius holes as singularities 4.3. Those
dot-sized holes are the same as the ones used in the computation of the
asymptotic energy of the drilled domains. The solution of this Neumann problem
is a single-valued scalar field that is strongly linked to the multivalued
angles of a cross field. Indeed, we can define an angle for each cross of a
smooth cross field. The computed single-valued field and this phase field have
perpendicular gradients and form together a complex harmonic field almost
everywhere. This Neumann problem could lead to a new simpler way to represent
cross fields and help to extract useful size information from it. Its inputs
consist of configurations of both the inner and the boundary singularities.
This last method offers more freedom to a user than the main method of our
paper while this last one creates a boundary singularity configuration that
minimizes our generalized Ginzburg-Landau energy for piecewise-smooth domain.

##  1 Related work

Let us briefly examine the literature that is either directly related to our
work or that follows the same objectives. For broader literature reviews see
the papers of Vaxman et al. [9] and Bommes et al. [10, 5]. Concerning the
theory on mesh generation in two dimensions, we refer the reader to a paper of
Bunin [11].

In the literature, there are two main ways to generate a block-structured quad
mesh. Either a cross field is computed to extract a quad mesh via a
parameterization or the parameterization is directly computed. However, there
is a strong links between the two approaches [12]. The equivalent of
singularities in the context of cross fields are cones in the context of
direct parameterization. In either case, the quality of the resulting object
depends on the amount of distortion in it, the number of singularities/cones,
and the alignment with features and boundaries [12]. The methods based on
cross fields always emphasizes the boundary alignment whereas parameterization
methods are often used for domain without boundaries [13] or without boundary
alignment (e.g. for texture mapping [14]). There are integer constraints in
parameterization to ensure the possibility to extract a block-structured
quadrangular mesh generation from it. These constraints are called
quantization and form a complex combinatorial problem (see [6] for a way to
solve them). Thus, the first quality that we seek in cross field for block-
structured quadrangular meshing is the alignment with the domain boundaries.
On one hand, only singularities that make sense in a quadrangular mesh are
used in the cross field context, i.e. singularities that are such that their
degree implies an integer number of neighbors for the corresponding node in
the underlying mesh. On the other hand, in parameterization, the correct
connection is not always requested and depends on the use of integer
parameters. In parameterization for four-sided zone construction, we have two
perpendicular fields and if we want an alignment to the domain boundary, an
iso value of one of the two fields must correspond to the boundary. With
regard to the cross fields, alignments with boundaries or features are
straightforward. Indeed it just imply that the tangents of the boundary or
feature curves corresponds to one of the vectors of the cross field at each
point. When we extract a parameterization from a cross field, the cones are
located at field singularities and gradients of parametric coordinates are
aligned with the field vectors [12].

Myles and Zorin [12] have proposed a method that place singularities to
minimize the size distortion in a parameterization. Regarding boundary-aligned
parameterization, they show that metric distortion may be reduced by cone
chains but that there will always have a trade-off, which can be arbitrated by
a regularization factor, between the distortion and the number of cones. In
our context, this paper therefore confirms that adding singularities can help
to limit the difference in size on computed field or parameterization and
therefore on block-structured quadrangular mesh extracted from them. The same
tradeoff appears in the Ginzburg-Landau theory. It can be controlled either by
a factor $\epsilon$ or by the radius of singular holes. It is important to
have a tradeoff because placing a lot of singularities can lead to the
appearance of very small blocks compared to the average block size in the
block-structured mesh that we can extract from a cross field. As these factors
tend to zero, we obtained asymptotic Ginzburg-Landau energy for which the main
minimization concerns the number and order of singularities. Contrary to what
we propose, the method Myles and Zorin [12] does not allow a user to impose
singularities. One point that seems particularly interesting to us is that
they include metric distortion directly in the optimization energy. The paper
of Myles and Zorin [12] does not consider additional integer constraints
needed for quadrangulation.

In addition to cross-field-based methods for constructing parameterization,
for the purpose of block-structured mesh generation, methods based on
streamlines are common. These streamlines extend the direction adjacent to
field singularities. Contrary to parameterization, streamline tracing methods
have a direct connection with the topology of the cross field. Leveraging the
Ginzburg-Landau theory, Viertel and Osting [8] claim that Ginzburg-Landau
cross field generation methods can be used to produce a cross field whose
separatrices divide the domain into four sided regions. They investigate the
mathematics of cross field generation and cross field guided quad meshing via
streamline tracing. They also made the observation that many of the
computational methods currently used for cross field design attempt to
minimize an energy close to the Ginzburg-Landau energy. They also highlight
the difficulty to compute numerically asymptotic minimizers of the energetic
Ginzburg-Landau functional. They suggest that this difficulty is link with the
fact that the Ginzburg-Landau energy assigns too much weight to singularities.
They further emphasize a slight influence of the refinement of the initial
discretization in the discretely computed energy. An important point of this
paper is that Viertel and Osting use ad-hoc methods to define boundary
singularities and hope to see the Ginzburg-Landau theory extended to domains
with piecewise-smooth boundaries. Precisely, later in this paper, we propose
such an extension of the Ginzburg-Landau theory.

Beaufort et al. [7] propose to directly solve the energetic Ginzburg-Landau
functional with a finite element formulation and with the use of edge-based
Crouzeix-Raviart interpolation functions. They also give theoretical results.
Particularly, they show that the positions of the singularities in the
computed cross fields coincide with the positions of irregular vertices in the
mesh that can be extracted from it.

Ray et al. [3] give a definition of the singularities of N-symmetry direction
fields close to the one that we propose on this article. The proposed
definition allows to relate singularities to the topology of the surface.
Specifically, they provide a generalization of the Poincaré-Hopf theorem to
N-symmetry direction fields on 2-manifolds. The proposed algorithm allows to
produce N-symmetry directional fields based on user-defined sets of
singularities. The approach of Ray et al. uses a greedy algorithm to constrain
singularities through the topology. The variables used to represent the
directions are angles relatively measured to a given arbitrary direction. This
implies an inherent difficulty since the fields of phases are often
multivalued. The researchers formulate the generation of cross fields as a
quadratic form minimization solved via conjugate gradient. However, no quality
estimations of the user-based locations of singularities are given. A lot of
topological analysis of interest are done in the article of Ray et al. [3]. In
direct link with the work that we will present, Ray et al. introduce the idea
of replacing singularities with holes, so that they can be handled as borders.
However, this paper present a different approach, as our formulation has
strong links with Ginzburg-Landau theory. Furthermore, our work is not based
on angles, which imply harder numerical computations. Last but not least, we
work with an efficient linear problem with some very localized quadratic
constraints.

Crane et al. [15] have developed a particular efficient way to impose
singularities in direction fields. They based their method on consistent
connection and the dual mesh of the input triangular mesh. Such triangular
meshes are classically used as a base in direction field generation. We also
use them as a base for all the numerical computations that we described in our
article. Crane et al. demonstrate that such consistent connections can be
computed easily and efficiently. They use homology consistent constraints on
particular cycles, contractible ones and non-contractible ones, to construct
their consistent connections. Such as for the cones in parameterization, they
concentrate curvature on singularities but in the same time keep a consistent
parallel transportation law. Their method summarizes in just finding the
minimum 2-norm solution of a constrained linear system. Their utilization of
connection on dual meshes gives interesting links between smooth and discrete
geometry. Just like us, they highlight the importance of topological
constraints. An important point is the enlightened tradeoff between total
control over the singularities and simplicity of the methods for singularity
placements. A remarkable point is that their method allows to never generate
non-wanted singularities. However, their method does not focus on boundaries
or features alignment that are essential in the purpose of block-structured
quadrangular mesh generation for high performance numerical simulation.

##  2 Theoretical background

A cross field is a particular case of a $n$-symmetry directional field with
$n=4$ (for an historical paper on this matter see [3]). It can be seen as a
quadruple of same norm vectors forming a regular cross, i.e., that are either
orthogonal or opposite to each other. A cross field $c$ on a planar region $G$
is a mapping that associates a cross $c(P)$ to every point $P\in G$. The four
vectors composing a cross $c(P)$ can be regarded as the $4^{th}$ roots of a
unique complex number $u(P)=c^{4}(P)$ (for a justification see [7]). This
representative complex number (or 2-d vector) is in general not parallel to
any of the branches of the cross it represents. A smooth cross field can thus
be represented unambiguously as a continuous $u\,:G\mapsto\mathcal{C}$ complex
function. The cross field $c$ is then the field of the $4^{th}$ roots of $u$
as illustrated in the Figure 1.

Fig 1: The crosses of the cross field represented on the left image are the
$4^{th}$ roots of the vectors of the vector field represented on the right
image.

We will take advantage of this link to build discrete cross field analysis on
the very established and rich continuous vector field analysis. In particular,
one singularity of the representing vector field will correspond to four
singularities of the derived cross field. These four singularities will have
an order equals to one fourth of the one of the vector field (for a
justification see [7]).

Vector fields and complex functions are equivalent concepts on planar regions
and we shall use indifferently one or the other according to the context. We
will see that the notions of degrees and indices are strongly linked and that
they allow, via the Poincaré-Hopf theorem, to characterize the total index of
a field on the basis of the Euler-Poincaré characteristic of the domain in
which it is calculated.

###  2.1 Brouwer degree

The Brouwer degree $\text{deg}(v,\gamma)$ of a vector field $v$ on a smooth
closed curve $\gamma$ plays a central role in the context of this paper. In
mathematical terms, it is defined as the degree of the application $v$
regarded as a map from $\gamma$ into $S^{1}$, where $S^{1}$ represents the
unit circle in the complex plane. Such an application is not always defined
which will in those cases be the same for the Brouwer degree. Nevertheless,
said application and Brouwer degree are well-defined whenever the field $v$
and the curve $\gamma$ are locally planar, so that a canonical mapping
$S^{1}\mapsto S^{1}$ can be found. In those instances, the Brouwer degree has
a simple and practical geometrical interpretation as the value of the integral

$\text{deg}(v,\gamma)=\frac{1}{2\pi}\int_{\gamma}d\theta(v)$ (1)

where $\theta(v)$ is the angle formed by the vector field $v$ and a fixed
reference vector, e.g., the x-axis of a local planar coordinate system. Note
that equivalent formula without the definition of an angle exist (for
definitions either via the determinant of a matrix, via an exterior product or
via lifting see [16]). That angle is counted positively according to the
orientation of the curve $\gamma$. The Brouwer degree is integer-valued for
vector fields and it is independent of the orientation of $\gamma$.

###  2.2 Index of a singularity

${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}\gamma}{\color[rgb]{0.82,0.01,0.11}{}_{1}}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}\gamma}{\color[rgb]{0.82,0.01,0.11}{}_{2}}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}\gamma}{\color[rgb]{0.82,0.01,0.11}{}_{3}}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}\gamma}{\color[rgb]{0.82,0.01,0.11}{}_{4}}$
Fig 2: The first three vector fields represented in black have a Brouwer
degree equals to 1 around the red paths. We have a source on $\gamma_{1}$, a
sink on $\gamma_{2}$ and a vortex on $\gamma_{3}$. The Brouwer degree of the
fourth vector field $v$ is -1 on $\gamma_{4}$ (saddle).

A first situtation where the Brouwer degree exists is when the curve $\gamma$,
defined on a smooth surface $M$ ($\gamma\subset M$), is infinitesimal, and
therefore can be considered locally planar. That implies that the presented
method could easily be extended to manifolds. Let $v$ be a vector field
defined on $M$. The index of an isolated singularity of $v$ located at a point
$P_{i}\in M$, is defined as the Brouwer degree

$\text{index}(P_{i})=\text{deg}(v,\gamma_{i}(P_{i}))=\frac{1}{2\pi}\int_{\gamma_{i}(P_{i})}d\theta(v),$
(2)

where $\gamma_{i}\subset M$ is a circle centered at $P_{i}$ of infinitesimal
radius, so that it encloses no other singularity. The index is 1 if the
singularity is of the type source, vortex or sink, and it is -1 for a saddle
type singularity (see figure 2).

###  2.3 Poincaré-Hopf theorem

The notion of index entails a very fruitful characterization for the isolated
singularities of vector fields (and also for cross fields as we shall see
further). However it has also a topological significance culminating in the
Poincaré-Hopf theorem (see [3, 7] for adaptations of this theorem in the
context of cross field). This theorem states that for a vector field $v$ with
isolated singularities defined on a closed differentiable manifold $M$ one has

$\sum\limits_{i=1}^{N}\text{index}(P_{i})=\chi(M)$ (3)

where $\chi(M)$ is the topological Euler-Poincaré characteristic of the
manifold $M$. In dimension 2, the characteristic of $M$ is given by

$\chi(M)=2-2g-b$ (4)

where $b$ is the number of connected components of the boundary $\partial M$,
and $g$ is the genus of the surface, i.e., the maximum number of cuttings
along non-intersecting closed curves that will not make the surface
disconnected.

###  2.4 The topological constraint for planar vector fields on smooth domains

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\Gamma}{\color[rgb]{0,0,0}{}_{ext}}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}G}$$\Gamma_{2}$$\Gamma_{1}$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\Gamma}{\color[rgb]{0,0,0}{}_{ext}}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}G}{\color[rgb]{0.82,0.01,0.11}{}_{\rho}}$$\Gamma_{2}$$\Gamma_{1}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}\gamma}{\color[rgb]{0.82,0.01,0.11}{}_{1}}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}\gamma}{\color[rgb]{0.82,0.01,0.11}{}_{2}}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}\gamma}{\color[rgb]{0.82,0.01,0.11}{}_{3}}$

Fig 3: Illustration of $G$, $\Gamma_{ext}$, the $\Gamma_{i}$ and the
$\gamma_{j}$.

The Poincaré-Hopf theorem will be developed so that only degrees and indices
will be paired, taking advantage of the fact that we work in planar. This will
lead to the definition of a global constraint that we will use to construct
our method for placing singularities within a domain.

The Poincaré-Hopf theorem highlights a profound link existing between a purely
analytical concept (the zeros of a vector field) and a purely topological one
(the Euler-Poincaré characteristic of the domain of definition of that vector
field). This link is extremely important in practice as it restricts the set
of configurations (degrees and positions) of singularities that a surface $M$
can contain based on its Euler-Poincaré characteristic $\chi(M)$. The theorem
is however only applicable to closed surfaces, or, if the surface has a
boundary with fields $v$ normal to that boundary. For the purposes of this
paper, a new interpretation of the generalizations of (3), valid for wedge-
shaped geometries, is developed later 4.2.

Let $G\subset\mathbb{R}^{2}$ be a planar region, whose boundary

$\partial G=\Gamma_{\text{ext}}-\Gamma_{1}-\cdots-\Gamma_{N_{\text{int}}}$

is composed of one external boundary $\Gamma_{\text{ext}}$ oriented counter-
clockwise, and of $N_{\text{int}}=b-1$ internal boundaries $\Gamma_{k}$
oriented clockwise (see Figure 3). The genus of a planar surface being zero,
(3) writes now

$\sum\limits_{i=1}^{N}\text{index}(P_{i})=2-b=1-N_{\text{int}}$ (5)

for planar surfaces. On these surfaces, the Brouwer degree of a vector field
$v$ can be defined for any closed curve $\gamma$, even non-infinitesimal.
Indeed, a parametrisation of the curve and a global coordinate system then
allows representing the vector field on any curve as a $S^{1}\mapsto S^{1}$
mapping. The Brouwer degree $\text{deg}(v,\gamma)$ is zero if $v$ is smooth in
the interior of $\gamma$, and it is otherwise equal to the sum of the indices
of the singularities enclosed in $\gamma$.

For smooth boundaries, if the vector field $v$ is normal to $\partial G$ one
can then evaluate the degrees $\text{deg}(v,\Gamma_{\text{ext}})=1$ and
$\text{deg}(v,\Gamma_{k})=-1$ for $k=1,\dots,N_{\text{int}}$, so that the
relationship

$\text{deg}(v,\Gamma_{\text{ext}})=\sum\limits_{i=1}^{N}\text{index}(P_{i})+\sum_{k=1}^{N_{\text{int}}}\text{deg}(v,\Gamma_{k})$
(6)

holds as an alternative expression of the Poincaré-Hopf theorem for planar
vector fields (see Figure 3).

It has been proven by Bethuel et al. [17] that this relationship remains valid
if the field $v$ has arbitrary degrees $\text{deg}(v,\Gamma_{\text{ext}})$ and
$\text{deg}(v,\Gamma_{k})$ on the boundary of $G$.

###  2.5 Planar cross field and Ginzburg-Landau functional

Let us present the mathematical theory of Ginzburg-Landau adapted to the cross
fields.

We will consider any Euler-Poincaré characteristic and require that one branch
of each boundary cross is parallel to the local outer normal $\mathbf{n}$ as
boundary conditions. Note that this condition shall not be respected on
corners of piecewise smooth boundary but the implication in terms of degrees
will be tackled via the definition of boundary singularities on such corners.
For regular boundary crosses (that are present almost everywhere on our
piecewise smooth domains of interest), the boundary condition writes

$c(P)\cdot\mathbf{n}\in\\{-1,0,1\\},\forall
P\in\\{\Gamma_{\text{ext}},\Gamma_{1},\ldots,\Gamma_{N_{\text{int}}}\\},$ (7)

i.e., two branches of a unit boundary cross are aligned with the boundary and
the other two branches of this unit boundary cross are opposite and parallel
to the local unit outward normal to the boundary.

Let’s temporary assume that we are working on smooth domains. If
$g_{c}\,:\\{\Gamma_{\text{ext}},\Gamma_{1},\ldots,\Gamma_{N_{\text{int}}}\\}\rightarrow
S^{1}$ is the boundary condition, smooth on each boundary $\Gamma_{i}$, then,
one observes that its Brouwer degree is fixed,
$\deg(g_{c},\\{\Gamma_{\text{ext}},\Gamma_{1},\ldots,\Gamma_{N_{\text{int}}}\\})=1-N_{\text{int}}$.

However, the branches of a cross are redundant as it is enough to know any one
of them to reconstruct the whole cross. At the same time, the representation
of a cross by one of its branches is also ambiguous as there exist no
systematic rule to select the representative branch of a cross. For these
reasons, it is convenient to solve indirectly for the field $u=c^{4}$. The
$4^{th}$ power maps all four branches of the cross onto a unique complex
number, which thus offers a smooth and unambiguous representation for the
cross. The boundary condition becomes for this vector field
$g=g_{c}^{4}\,:\Gamma_{\text{ext}}\rightarrow S^{1}$, with the Brouwer degree

$\deg(g,\\{\Gamma_{\text{ext}},\Gamma_{1},\ldots,\Gamma_{N_{\text{int}}}\\})=4-4N_{\text{int}}.$
(8)

Consequently, to compute a cross field aligned with the domain boundaries and
minimizing the Ginzburg-Landau functional, one is led to solve the variational
problem

$\min\limits_{H^{1}_{g}(G,\mathbb{R}^{2})}\int_{G}|\nabla
u|^{2}+\frac{1}{4\epsilon^{2}}(|u|^{2}-1)^{2}.$ (9)

with $H_{g}^{1}=\\{u\in H^{1}(G,\mathbb{R}^{2});\;u=g\text{ on
}\\{\Gamma_{\text{ext}},\Gamma_{1},\ldots,\Gamma_{N_{\text{int}}}\\}\\}$. The
gradient term ensures smoothness, while the penalty term ensures $|u|$ to be
as close to 1 as possible. The minimized term of (9) is called the Ginzburg-
Landau energy of the vector field $u$. If $\epsilon$ is small, the minimizers
$u_{*}$ will therefore be close to 1 outside a small set of isolated
singularities, whose cumulated indices sum up to $4-4N_{\text{int}}$.
Mathematically, when we minimize (9) as $\epsilon\rightarrow 0$ there is talk
of asymptotic minimizers of the energetic Ginzburg-Landau functional (9). In
the discrete model, singularities are regions of size $\epsilon$ where $|u|$
deviates significantly from 1 towards zero. Note that an adaptation of this
energy for two-dimensional manifolds that only implies linear added terms
(linked to the Christoffel symbols) has been proposed in [18].

Now (6) and (8) imply that

$\sum\limits_{i=1}^{N}\text{index}(P_{i})=\text{deg}(v,\\{\Gamma_{\text{ext}},\Gamma_{1},\ldots,\Gamma_{N_{\text{int}}}\\})=4-4N_{\text{int}}.$
(10)

If we count in the points $P_{i}$ the internal singularities and also the
boundary singularities, the formula (10) holds for domains with piecewise
smooth boundary.

Functionals of the form (9) were originally introduced by Ginzburg and Landau
in the study of the phase transition problems occurring in superconductivity
(see an overview of their research in [19]).

The approach of Beaufort et al. (see [7]) sees the optimal singularity
configurations as the zeros and associated degrees of the minimizers of
$E_{\epsilon}$ for $\epsilon$ small enough. Unfortunately, this formulation
can not numerically ensure that only a minimum number of singularities,
necessary to fulfill (3), appears and thus is not always able to reproduce the
asymptotic mathematical results obtained in [17].

Fig 4: Here are some results obtained via the method described in [7] based on
the Ginzburg-Landau functional (9) for rings with different ratios of the
small radius $r_{1}$ to the large radius $r_{2}$. For a parameter
$\epsilon=1/100$, we obtain, on the left figure (ring with $r_{1}/r_{2}=0.4$),
8 singularities (four of degree equals to 1/4 and four of degree equals to
-1/4) in our cross field and on the right figure (ring with
$r_{1}/r_{2}=0.46$) zero singularity. For a given $\epsilon$, when the ring
becomes too thick, singularities appear.

The ring is an example where this phenomenon can occur (see figure 4).
Although its Euler-Poincaré characteristic is $0$, numerically, singularities
appear if the ratio of the small radius to the great radius is lower or equal
to $0.4$ for $\epsilon=1/100$ (see figure 4). With a smaller $\epsilon$ we may
have converged to the radial solution with no singularities. However, the use
of $\epsilon$ is a problem in itself: the $\epsilon$ needed to obtain a
numerical result analogue to the asymptotic mathematical minimizer of (9)
depends strongly on the domain and must possibly be very small, which may be
hard to work with numerically.

On the other hand, the supplementary singularities can be viewed as an
advantage for the quadrangular mesh generation because they can be used to
limit the size distortions of an elementary quad cutting. Let us compare the
two rough quad components of Figure 5. For thick ring, the radial quad
decomposition of the right figure of 5 has great size difference between its
elements. On the contrary, the singularities of opposite sign in the left
figure of 5 enable to thwart those deformations.

Fig 5: Obtained meshes via rough directions propagation (priority is given to
propagating the directions adjacent to the singularities) for the cross fields
of figure 4. For a parameter $\epsilon=1/100$, we obtain, on the left figure
(ring with $r_{1}/r_{2}=0.4$), a quad mesh with less deformations (especially
if we refine it) than the quad mesh of the right figure (ring with
$r_{1}/r_{2}=0.46$) that contains no singularity.

Judging which of the two cross fields is best from a numerical point of view
depends on the properties sought after - for example, emphasis on preferential
directions, or, on the contrary, search for uniformity of the elements.
Furthermore, the addition of singularities should potentially make it possible
to design mesh size for the block-structured quadrangular meshes resulting
from the computed cross fields.

Here we have a working but rigid method to construct cross fields on surfaces.
In the section 3 we will develop an alternative method in order to choose the
configurations of inner singularities.

###  2.6 How to evaluate the choice of a singularity location

To find the optimal singularity configuration for various radius of the
singular holes with the method presented in 3, we propose the following energy
that is strongly linked with the Ginzburg-Landau energy. This energy is the
energy of the scalar field that is solution of the following problem:

$\begin{cases}-\Delta\Phi=0&\text{ in }G_{\rho},\\\
\Phi=\text{Const.}=C_{i}&\text{ on }\gamma_{i},\;i=1,2,\ldots N,\\\
\Phi=0&\text{ on }\partial G,\\\
\int_{\gamma_{i}}\frac{\partial\Phi}{\partial\nu}=2\pi
d_{i}&i=1,2,\ldots,N,\end{cases}$ (11)

where $G_{\rho}$ is the domain $G$ in which we have drilled holes of radius
$\rho$ with boundaries $\gamma_{i}$ such as $\chi(G_{\rho})=0$ and where $\nu$
is the outward normal to the $\gamma_{i}$ (see Figure 3 for an illustration).
Here the $C_{i}$ are not given but are unknown constants that are part of the
problem (see [17]).

Theorem I.1 in [17] tells us that

$E_{\rho}=\int_{G_{\rho}}|\nabla\Phi|^{2}=\inf\limits_{v\in
H^{1}(\mathbb{R}^{2},G_{\rho})}\int_{G_{\rho}}|\nabla v|^{2}.$

That means that the solution of the proposed problem have the same energy as
the one of the regular cross field proposed in 3.

As the problem (11) is a simple linear Neumann problem, it gives a simple way
to cheaply compares singularities configurations that could be used with the
main method of this paper that will be presented in the beginning of the
section 3.

But, if we are not interested in the impact of the singularity radius or if we
let those radius tend to zero, we could even compare the singularity in term
of renormalized Ginzburg-Landau energy by direct computation.

####  2.6.1 The unifying renormalized Ginzburg-Landau energy

As we take the $\gamma_{i}$ as being disk of radius $\rho$, we have

$\lim\limits_{\epsilon\rightarrow
0}\min\limits_{H^{1}_{g}(G,\mathbb{R}^{2})}\frac{1}{2}\int_{G}|\nabla
v|^{2}+\frac{1}{4\epsilon^{2}}\int_{G}(|v|^{2}-1)^{2}=\lim\limits_{\rho\rightarrow
0}\inf\limits_{v_{\rho}\in\mathcal{E}}\int_{G_{\rho}}|\nabla u|^{2}=+\infty$

where

$\mathcal{E}=\left\\{v\in H^{1}(G_{\rho};S^{1})\left|\begin{aligned} v=g\text{
on }\partial G\text{ and }\\\ \deg(v,\gamma_{i})=1\;\;\forall
i\end{aligned}\right.\right\\}.$

But as explained in [17], there is a way to compare energetically the two
approaches corresponding with those energies by removing the infinite core
energy from them. The resulting common energy is called renormalized energy.

Theorem I.2 and Theorem I.7 of [17] prove that there exists a unique minimizer
$u_{\rho}$ for the problem

$\min\limits_{u\in\mathcal{E}}\int_{G_{\rho}}|\nabla u_{\rho}|^{2}$ (12)

and that the following expansion holds:

$\frac{1}{2}\int_{G_{\rho}}|\nabla u_{\rho}|^{2}=\pi
d|\log\rho|+W(\\{P\\})+O(\rho)\text{ as }\rho\rightarrow 0$

where $W$, called the renormalized energy, stays bounded as $\rho\rightarrow
0$ and where $\\{P\\}=(P_{1},P_{2},\ldots P_{d})$ is the configuration of the
centers of the circles $\gamma_{i}$. We call singular “core energy” $\pi
d|\log\rho|$. Other examples of removing singular core energy appear in
physics (see [20] for another example).

The renormalized energy $W$ does not depend on $\rho$ and is the same as the
renormalized energy of the energetic Ginzburg-Landau functional. More
specifically, it is possible to computed the common renormalized energy $W$
via a set of point locations and degrees rather than via a function such as in
(9). Indeed the positions and degrees of the singularities of the minimizing
field of the energetic Ginzburg-Landau functional, which is defined for maps
$u\in H^{1}(G;\mathbb{C})$,

$E_{\epsilon}(u)=\frac{1}{2}\int_{G}|\nabla
u|^{2}+\frac{1}{4\epsilon^{2}}\int_{G}(|u|^{2}-1)^{2}$

or the positions and degrees of the singularities that minimize

$\min\limits_{u\in\mathcal{E}}\int_{G_{\rho}}|\nabla u_{\rho}|^{2}$

are the same as the positions of the points of the configuration
$b=(b_{1},b_{2},\ldots,b_{d})$ that minimize the renormalized Ginzburg-Landau
energy, as $\epsilon\rightarrow 0$ and if $d=\deg(g,\partial G)$,

$W(b)=-\pi\sum\limits_{i\neq j}\log|b_{i}-b_{j}|+\frac{1}{2}\int_{\partial
G}\Phi(g\times g_{\tau})-\pi\sum\limits_{i=1}^{d}R(b_{i})$

where $\Phi$ is the solution of the linear Neumann problem

$\begin{cases}\Delta\Phi=2\pi\sum\limits_{i=1}^{d}\delta_{b_{i}}&\text{ in
}G,\\\ \frac{\partial\Phi}{\partial\nu}=g\times g_{\tau}&\text{ on }\partial
G,\end{cases}$

where $\nu$ is the outward normal to $\partial G$ and $\tau$ is a unit tangent
vector to $\partial G$ such that $(\nu,\tau)$ is direct, the $\delta_{b_{i}}$
corresponds to the singularity degrees and

$R(x)=\Phi(x)-\sum\limits_{i=1}^{d}\log|x-b_{i}|.$

This renormalized energy teaches us that whether a formulations uses holes
with a certain radius or the factor epsilon, they all can be compared in terms
of the asymptotic Ginzburg-Landau energy by providing their singularities
positions and degrees. Indeed, those results show us that the minimal
configuration of holed singularities from (13) as $\rho\rightarrow 0$ is the
same as the configuration of the zeros of the minimizers of (9). Furthermore,
they give a way to compare different choices of singularity configuration in
terms of renormalized Ginzburg-Landau energy. In particular, Theorem 0.2 of
[17] proves that the zeros of minimizers of (9) converge to minimizers of $W$.

A scalar field adapted from $\Phi(x)$ to work with cross field is solved in
the section 4.3 and can be used to compute the proposed renormalized energy.
One element that makes this a particularly simple problem is that the
equivalent of our holes materializes as dots and the equivalent of our degrees
as deltas.

Without numerical computation, the analysis of the renormalized energy $W$
teaches us the following things that could be use to compare singularity
configurations. First, $W\rightarrow+\infty$ as two singularities $P_{i}$,
$P_{j}$ of both positive or negative degrees coalesce Secondly,
$W\rightarrow+\infty$ as a $P_{i}$ tends to $\partial G$. Third, $W$ decrease
as one singularity of positive degree and one singularity of negative degree
coalesce. Those properties imply that numerically singularities of
configuration minimizing the Ginzburg-Landau energies (the renormalized ones
and the one of (9)) will rather take position moderately far to concave
boundaries, more far from convex boundaries and will repulse or attract each
other following their degree signs.

While the classical energetic Ginzburg-Landau functional involves finding a
function representing a vector field that minimizes a functional, the second
consists only in finding a configuration of points as a minimizer. Moreover,
the calculation of $\Phi$ in the renormalized energy can be done for any
configuration of singularities even though it would not minimize the Ginzburg-
Landau renormalized energy. It is therefore a major lead for efficient cross
field computations.

###  2.7 Pinning strategy

Let us now consider another way of looking at the imposition of singularities.
This kind of method is of great interest in physics but is less interesting in
the context of cross fields, as we will explain briefly.

The pinning strategy consists in modifying the Ginzburg-Landau functional in
order to pin singularities in desired places. We can either change the non-
linear term of this functional or add a weight that varies on the domain in
front of the gradient term. In physics this type of approach is called
pinning. For more information on changes on the penalty term see [21] and [22]
and on weight additions see [23].

Numerically, the problem is that those methods are even hard to solve than the
method [7].

For cross fields, a classical formulation for their representation vector
could be the following. We are looking for $u\in H^{1}(G,\mathbb{R}^{2})$ that
minimizes the functional

$\int_{G}|\nabla
u|^{2}+\frac{1}{4\epsilon^{2}}\int_{G}\left(f(x)-|u|^{2}\right)^{2}$

with $\deg(u,\partial G)=4$ and with

$f(x)=\exp\left(\|(\sqrt{x^{2}+y^{2}}-0.5)\|\right)$

for a pin on a circle centered at the origin of radius $0.5$.

###  2.8 Quad mesh generation from cross fields

There are a lot of different methods for extracting quadrangular meshes from
cross fields. Among the numerous utilizations, we can highlight 3 interesting
approaches. The first one consists in generating streamlines from
singularities and to use the underlying cutting to compute quad mesh (see, as
examples, [24] and [25]). The second works with parametrization via scalar
fields and mixed integer solvers (see, as examples, [26] and [27]). The third
approach uses Morse parameterization hybridization (see [28]).

##  3 Materials and methods

We will begin by presenting our method for placing singularities in a domain.
We will see that this method does not allow us to explicitly set boundary
singularities. The singularities that will appear on the boundaries will in
fact be the ones that minimize a new energy that we construct in section 3.2.
This will imply that it is not easy to know how to respect the Poincaré-Hopf
theorem since this theorem takes into consideration both inner and boundary
singularity degrees. For the practical use of our method, we will propose a
new interpretation of the generalizations of the Poincare-Hopf theorem in
section 4.2.

###  3.1 Singularity placements via perforation strategy

As we have seen, regular vector fields can be computed on any surface as long
as we have drill holes with prescribed degrees on their boundaries such as to
respect the Poincaré-Hopf theorem. Furthermore, the analysis of Ginzburg-
Landau energy have lead to simplified formulations that allows to have non-
minimal number and non-minimizing (in terms of Ginzburg-Landau energy)
configuration of singularities (see [17]).

####  3.1.1 Description of the problem

The first idea developed in this paper, following a similar approach presented
in [17], is that the unavoidable singularities could be trapped in tiny holes,
i.e., excluded from the computational domain. In order to do so, suitable
conditions must be prescribed on the boundaries of these holes. One would then
be left with the computation of a regular (singular free) field that is
numerically cheap and stable.

Let thus $B(P_{i},\rho)$ be $N$ small holes drilled in the domain $G$,
$\gamma_{i}=\partial B(P_{i},\rho)$, and
$G_{\rho}=G\backslash\bigcup_{i}B(P_{i},\rho)$ be the computational domain
excluding the holes. The points $P_{1},P_{2},\ldots,P_{N}$ at which the holes
are centered can be placed freely.

In [17] there is a proof of the existence of minimizers $v_{\rho}$ for the
problem

$\min\limits_{v\in\mathcal{E}}\int_{G_{\rho}}|\nabla u|^{2}.$ (13)

where

$\mathcal{E}=\left\\{v\in H^{1}(G_{\rho};S^{1})\left|\begin{aligned} v=g\text{
on }\partial G\text{ and }\\\ \deg(v,\gamma_{i})=1\;\;\forall
i\end{aligned}\right.\right\\}.$

The formulation we will use to place singularities will be slightly different
than (13). For the sake of clarity, the differences are detailed in the
following bulleted list.

  * •

The function $g$ will be such that $u$ has two of its $4^{th}$ roots parallel
to $\partial G$ so that the associated cross fields are tangent with $\partial
G$.

  * •

The used vector field will belong to $\mathbb{R}^{2}$ in the domain rather
than $S^{1}$ for cheap implementation. However, the boundary representation
vectors must not divert too much from $S^{1}$ in order to preserve a well-
defined notion of degree. If this notion is not well enough preserved, the
computation must be performed with an added constraint on their norms.

  * •

In the same way as for the previous bullet, we will suppose that $u\in S^{1}$
on the circles $\gamma_{i}$ and use the same notion of degree.

  * •

We will work on connected piecewise smooth bounded domain $G_{\rho}$.

  * •

The configurations of $\gamma=(\gamma_{1},\gamma_{2},...,\gamma_{p})$ will be
such as $d=\sum_{i=1}^{p}\deg(u,\gamma_{i})$ but without imposing that $p=d$
with $d=\deg(g,\partial G)$ for smooth domains. Furthermore, $d$ will be equal
to a slight generalization of the Brouwer degree for piecewise smooth domains
that we will give later on (see section 4.2).

Note that we will now have holes in our fields instead of the singularities of
the fields composed of unitary elements but we will keep the name
“singularity” since these holes always mimic classical singularities and
correspond also to singular nodes (nodes with strictly less or more than 4
neighbors) in the meshes derived from these fields.

Three advantages associated with this approach are worth mentioning. Firstly,
unlike what is shown in [7] and [8], a parameter $\epsilon$, whose value is
not clear and strongly dependent on the problem considered, is not needed
anymore. Note that the $\rho$ parameter has taken the place of the $\epsilon$
parameter but our problem stays more stable as long as we do not impose $\rho$
to be very small. Secondly, the use of holes to replace zero norm crosses as
singularities in $G$ implies that the computed field may be non-singular (or
regular) and smooth. The third advantage is the ability to better adapt to the
user’s needs and demands for quadrangular meshing.

However, one drawback of this approach is that, numerically, the work is
performed on reference vectors belonging to $\mathbb{R}^{2}$, while the
associated mathematical results have been obtained for vectors belonging to
$S^{1}$. We will then have a disparity on the norms of the crosses in
compensation for the rotations induced in the computed cross fields.

####  3.1.2 Discrete formulation and implementation

Numerically, calculations are performed in the domain via Crouzeix-Raviart
finite elements (see [29] for more details) on a triangular mesh generated via
GMSH (see [30] for more details). It is then a question of building a finite
element system. Its matrix of stiffness corresponds to a Laplacian within the
domain. At the same time, on each node at the boundary, we construct a cross
having two directions tangent to one of the two edges adjacent to the node.
Note that the implementation was done in Python using the GMSH python API
which can generate the Jacobians and shape functions of the initial triangular
mesh.

####  3.1.3 Definition of the constraint for holed singularities

Let $u\in H^{1}_{g}(G_{\rho},\mathbb{R}^{2})$ where $G$ is a piecewise smooth
bounded domain and $G_{\rho}$ equals $G$ perforated by $N$ circular holes of
boundaries $\gamma_{j},1\leq j\leq N$. The following constraint must be
respected

$\sum\limits_{j=1}^{N}\deg(u,\gamma_{j})=\chi\deg(g,\partial G)$ (14)

which is a particular form of the Poincaré-Hopf theorem.

We start from the definition of degree (1). It is important to note that, on
each $\gamma_{j}$, the variation of the angle could be defined as the
derivative of an arc-tangent. We have

$\displaystyle d\phi$
$\displaystyle=\left(\frac{1}{1+\left(\frac{v_{1}}{v_{2}}\right)^{2}}\right)\left(\frac{dv_{1}}{v_{2}}-\frac{dv_{2}v_{1}}{v_{2}^{2}}\right)$
$\displaystyle=\frac{v_{2}^{2}}{v_{1}^{2}+v_{2}^{2}}\left(\frac{dv_{1}}{v_{2}}-\frac{v_{1}}{v_{2}^{2}}dv_{2}\right).$

On $\gamma_{j}$, supposing that the elements of the field have norm 1 on each
circle $\gamma_{j}$, we have

$d\phi=v_{2}dv_{1}-v_{1}dv_{2}.$

However, the degree of $u$ on each $\gamma_{j}$ precisely equals to the closed
integral of $d\phi$ on $\gamma_{j}$ divided by $2\pi$. Thus, our general
constraint (14) can be reformulated as follows:

$\frac{1}{2\pi}\sum_{j=1}^{N}\int_{\gamma_{j}}d\phi=\chi\deg(g,\partial G).$

A formulation of the constraint imposed on a single hole will now be derived.
On each $\gamma_{j}$ we have local polar coordinates (see figure 6).

$r=\sqrt{v_{1}^{2}+v_{2}^{2}}\text{ and
}\phi=\arctan\left(\frac{v_{1}}{v_{2}}\right)$

where $v_{1}$ and $v_{2}$ are the Cartesian coordinates of the vector $v$
considered on $\gamma_{j}$.

$\theta$$\phi$$v$${\color[rgb]{0.25,0.46,0.02}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.46,0.02}\gamma_{j}}$
Fig 6: Definition of local polar coordinates.

With these new coordinates, $dr=0$ on $\gamma_{j}$. Furthermore, the
derivatives of our vector components only involve partial derivatives with
respect to the $\theta$ variable. For each $\gamma_{j}$, this leads to the
following equation:

$d\phi=\left(v_{2}\frac{\partial v_{1}}{\partial\theta}-v_{1}\frac{\partial
v_{2}}{\partial\theta}\right)d\theta,$

where $\theta$ is the angle to the center of each $\gamma_{j}$.

At this point, it is important to note that numerically, it will be necessary
to pay particular attention to the orientation on $\gamma_{j}$ in order to
truly match the general constraint (14).

In practical terms, the following centered differences will be worked with,
which can be understood by the figure 7:

$\displaystyle\left(\frac{\partial
v_{1}}{\partial\theta}\right)^{i}=\frac{v_{1}^{i+1}-v_{1}^{i-1}}{2\Delta\theta},$
$\displaystyle\left(\frac{\partial
v_{2}}{\partial\theta}\right)^{i}=\frac{v_{2}^{i+1}-v_{2}^{i-1}}{2\Delta\theta},$

$\displaystyle\left(\Delta\phi\right)^{i,i+1}$
$\displaystyle=\left(\frac{v_{2}^{i}\left(v_{1}^{i+1}-v_{1}^{i-1}\right)-v_{1}^{i}\left(v_{2}^{i+1}-v_{2}^{i-1}\right)}{2\Delta\theta}\right)\Delta\theta,$
$\displaystyle=\frac{1}{2}\left(v_{2}^{i}\left(v_{1}^{i+1}-v_{1}^{i-1}\right)-v_{1}^{i}\left(v_{2}^{i+1}-v_{2}^{i-1}\right)\right),$

where the indices placed in exponents indicate successive node numbers.

$\Delta\theta$$i+1$$i$$i-1$ Fig 7: Finite centered differences for
$\left(\frac{\partial v_{1}}{\partial\theta}\right)^{i}$ and
$\left(\frac{\partial v_{2}}{\partial\theta}\right)^{i}$.

Hence, for a hole $\gamma_{j}$ approximated by a regular polygon with $m$
edges, we have

$\sum\limits_{i=1}^{m-1}\left(\Delta\phi\right)^{i,i+1}\approx\int_{\gamma_{j}}d\phi.$

Finally, we have as a discrete approximation of the constraint equation (14):

$\frac{1}{2\pi}\sum\limits_{j=1}^{N}\left(\sum\limits_{i=1}^{m-1}\left(\Delta\phi\right)^{i,i+1}\right)_{\gamma_{j}}=\\\
\frac{1}{4\pi}\sum\limits_{j=1}^{N}\left(\sum\limits_{i=1}^{m-1}v_{2}^{i}\left(v_{1}^{i+1}-v_{1}^{i-1}\right)-v_{1}^{i}\left(v_{2}^{i+1}-v_{2}^{i-1}\right)\right)_{\gamma_{j}}\\\
\approx\chi\deg(g,\partial G).$

In practice, we will separate the overall condition into one condition on each
perforated hole in order to clearly choose the placement of the singularities.
The resolution of a $\Delta u=0$ with the imposition of conditions on the
boundary crosses to align with $\partial G$ leads to a linear system

$Kx=b$

with $x=(v_{1}^{1},v_{2}^{1},\cdots v_{1}^{n},v_{2}^{n})$.

We can rewrite the non-linear quadratic constraint as followed:

$x^{T}Mx=\sum\limits_{i=1}^{2n}\sum\limits_{j=1}^{2n}m_{ij}x_{i}x_{j}=2\pi\chi\deg(g,\partial
G).$

For example, if we have a single hole with nodes 4-38-39-40-41 on it, $M$ will
have the following shape:

$v_{1}^{4}$ $v_{2}^{4}$ $v_{1}^{38}$ $v_{2}^{38}$ $v_{1}^{39}$ $v_{2}^{39}$
$v_{1}^{40}$ $v_{2}^{40}$ $v_{1}^{41}$ $v_{2}^{41}$ $v_{1}^{4}$ -1 1
$v_{2}^{4}$ 1 -1 $v_{1}^{38}$ 1 -1 $v_{2}^{38}$ -1 1 $v_{1}^{39}$ 1 -1
$v_{2}^{39}$ -1 1 $v_{1}^{40}$ 1 -1 $v_{2}^{40}$ -1 1 $v_{1}^{41}$ -1 1
$v_{2}^{41}$ 1 -1

We then have

$x^{T}Mx=\;v_{2}^{4}v_{1}^{38}-v_{2}^{4}v_{1}^{41}-v_{1}^{4}v_{2}^{38}+v_{1}^{4}v_{2}^{41}+v_{2}^{38}v_{1}^{39}\\\
-v_{2}^{38}v_{1}^{4}-v_{1}^{38}v_{2}^{39}+v_{1}^{38}v_{2}^{4}+v_{2}^{39}v_{1}^{40}\\\
-v_{2}^{39}v_{1}^{38}-v_{1}^{39}v_{2}^{40}+v_{1}^{39}v_{2}^{40}+v_{1}^{39}v_{2}^{38}\\\
+v_{2}^{40}v_{1}^{41}-v_{2}^{40}v_{1}^{39}-v_{1}^{40}v_{2}^{41}+v_{1}^{40}v_{2}^{39}\\\
+v_{2}^{41}v_{1}^{4}-v_{2}^{41}v_{1}^{40}-v_{1}^{41}v_{2}^{4}+v_{1}^{41}v_{2}^{40}.$

In order to solve this problem, a Newton-Raphson scheme is computed with the
utilization of a Lagrange multiplier for each of the constraint of each of our
perforated holes. We then have each constraint $C_{\gamma_{i}}$ on each
$\gamma_{i}$ and its derivatives which are described by

$C_{\gamma_{i}}:=\lambda_{i}\left(x^{T}Mx-2\pi\deg(u,\gamma_{i})\right),$
$\frac{\partial
C_{\gamma_{i}}}{\partial\lambda_{i}}=x^{T}Mx-2\pi\deg(u,\gamma_{i}),$
$\frac{\partial C_{\gamma_{i}}}{\partial x}=2\lambda_{i}Mx,$
$\frac{\partial^{2}C_{\gamma_{i}}}{\partial\lambda_{i}\partial x}=2Mx,$
$\frac{\partial^{2}C_{\gamma_{i}}}{\partial\lambda_{i}\partial\lambda_{i}}=0,$

and

$\frac{\partial^{2}C_{\gamma_{i}}}{\partial x\partial x}=2M.$

The gradient of our problem becomes

$\tilde{b}=[b+2\lambda_{i}Mx,\;x^{T}Mx-2\pi\deg(u,\gamma_{i})]^{T}$

and the Hessian becomes

$H=\begin{bmatrix}K+2M&2Mx\\\ (2Mx)^{T}&0\end{bmatrix}.$

We solve this via Newton-Raphson scheme until convergence and if the Hessian
is singular at a particular iteration, we set for this particular iteration :

$\frac{\partial^{2}C_{\gamma_{i}}}{\partial\lambda_{i}\partial\lambda_{i}}=1.$

###  3.2 New computations for placing boundary singularities

When cross fields are computed on piecewise smooth closed domains, in addition
to the internal singularities we have already discussed, boundary
singularities may appear on the boundary corners. The corner angles to be
considered correspond to the total angle formed around the boundary corner
inside the domain. Thus, an outgoing corner will correspond to an angle
between $0$ and $\pi$ radians while an incoming corner will correspond to an
angle between $\pi$ and $2\pi$ radians.

$\displaystyle\alpha$ Fig 8: Comparison of the presence of a boundary
singularity of degree $1/4$ in red on the apex of the $\alpha(\approx\pi/2)$
angle or the absence of singularity in purple.

A boundary singularity implies a rotation of the boundary cross compared to
the local outward normal at the domain boundaries. For example, if around a
boundary corner, the border cross rotation is $\frac{\pi}{2}$ compared to the
local outward normal, we say that a $1/4$ boundary singularity has been placed
on this corner. The presence or the absence of a singularity on such an angle
implies different rotations for the computed cross field (see Figure 8) since
the field inside the domain must be smooth.

The boundary corners are usually considered as singularity for angle multiple
of k$\pi/2$, considering that a singularity of index 0 is in fact a regular
cross. For example the four corners of a square correspond to four 1/4
boundary singularities thus the Poincaré-Hopf theorem is respected without
internal singularities. The conventional distribution proposed in the
literature corresponds to the one given on Table 1 with different tolerances
around the given angles.

angles (in radians) | $2\pi$ | $\frac{3\pi}{2}$ | $\pi$ | $\frac{\pi}{2}$ | 0
---|---|---|---|---|---
index of the singularity | $-\frac{2}{4}$ | $-\frac{1}{4}$ | $0$ | $\frac{1}{4}$ | one boundary $\frac{1}{4}$ and one internal $\frac{1}{4}$
Table 1: Usual singularities angles distribution.
${\color[rgb]{0.25,0.46,0.02}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.46,0.02}\alpha\approx\pi}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}\beta\
\approx\
}{\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}\frac{\pi}{2}}$${\color[rgb]{0.25,0.46,0.02}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.46,0.02}\alpha\approx\pi}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}\beta\
\approx\
}{\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}\frac{\pi}{2}}$
Fig 9:  On the above left diagram, a singularity $1/4$ has been placed around
an angle $\alpha$ close to $\pi$. On the above right diagram, no singularity
has been placed around the same angle $\alpha$. We see that for this angle
$\alpha$, placing no, i.d. a zero degree, boundary singularity on it implies
less rotations on the computed cross field. Conversely, on the below left
diagram, a singularity $1/4$ has been placed around an angle $\beta$ close to
$\frac{\pi}{2}$. On the below right diagram, no singularity has been placed
around the same angle $\beta$. We see that for this angle $\beta$, placing a
1/4 boundary singularity on it implies less rotations on the computed cross
field.

To deal with boundary singularities there are two possibilities: to place the
singularities as we do with the problem described in the section about further
thoughts 4.3 or to focus only on the internal singularities but to count the
boundary singularities via some boundary angular jumps as we will do in the
section where we give a new interpretation of the Poincaré-Hopf theorem 4.2.

The angles proposed in the literature to place boundary singularities only try
to limit the rotation of the cross field around these boundary corners as it
is illustrated on Figure 9. However, to respect the Poincaré-Hopf theorem
while limiting the rotation on a computed cross field, we believe that
internal singularities must be taken into account.

In this section, informal calculations are used to justify the range of
boundary angles that are considered as singularities as well as their orders
depending on the boundary angles and the computational geometry. We believe
that the proposed energy to be minimized is a natural extension of the
Ginzburg-Landau energy for angular domains. This energy let the internal
singularities affect the angle ranges on which a boundary singularity must be
placed in order to minimize this boundary Ginzburg-Landau energy. The
classical k$\pi/2$ angles as boundary singularities belongs to the expected
singularity degrees. As we will see, the proposed ranges depend on whether
there are too many or too few boundary singularities to match the Euler-
Poincaré characteristic of the domain. Intuitively, if we have not enough
positive boundary singularities to match the Euler-Poincaré characteristic, we
must place boundary positive singularities on larger range of angles and vice-
versa. Keep in mind that the sums of internal and boundary singularities are
strongly linked by the Poincaré-Hopf theorem (3).

We will see three scenarios appear. In terms of implementation, the idea is to
compute the energy implied by those three scenarios and to choose the ranges
of angle for boundary singularities that correspond to the scenario of less
energy.

Note that for quad mesh purpose, boundary singularity of degree 1/2 are of no
interest and must be replaced by one 1/4 boundary singularity and one 1/4
internal singularity close to the boundary one.

####  3.2.1 Ginzburg-Landau energy for piecewise-smooth closed boundaries

Given a domain $\Omega$ with a piecewise $C^{1}$ boundary $\partial\Omega$, we
consider a cross field, whose elements are in $SO(2)/C_{4}$, where $C_{4}$ is
the rotation group of the square (square is preserved by rotation of
$\\{\frac{j}{4}2\pi|j\in\\{0,...,3\\}\\}$).

We assume that the angles formed at the boundary of $\Omega$ are
$\alpha_{1},\ldots,\alpha_{m}$ and taken between $0$ and $2\pi$. We let
$k_{1},\ldots,k_{m}$ be the indices of the singularities inserted in those
angles.

Note that this may include singularities of index $0$ which are just a regular
crosses. To respect the Poincaré-Hopf theorem, certain singularities have to
be placed in the domain. We let $l_{1},\ldots l_{n}$ be their indices.

On one hand, the singular energy induced by each featured angle $\alpha_{i}$
of the boundary is given by

$\sum\limits_{i=1}^{m}\frac{(\pi-\alpha_{i}-k_{i}2\pi)^{2}}{2\alpha_{i}}.$

The idea is that the rotational energy needed to not place a singularity on
some angle is given by the linear variation induced on the cross. This
variation is squared and normed with the size of the angle in order to
construct an energy. The $1/2$ term is arbitrary and stems from the fact that
we described some fictive kinetic energy. This way, placing a 1/4 singularity
on a flat edge will cost: $\frac{\pi^{2}}{4}\frac{1}{2\pi}=\frac{\pi}{8}$ in
terms of our proposed Ginzburg-Landau like energy. This energy also implies
that put no singularity on a flat angle or a 1/4 singularity on a $\pi/2$
angle or a -1/4 singularity on a $3\pi/2$ angle of the external boundary will
result in 0 singular energy (see Figure 8).

On the other hand, the singular energy induced by all the inner singularity is
given by

$\frac{1}{2\pi}\left(\sum\limits_{i=1}^{n}\left(2\pi|l_{i}|\right)^{2}\right).$

since this is the total rotational energy, computed as above, created in the
domain. This way, placing an extra 1/4 singularity in the domain will cost:
$\frac{1}{2\pi}\left(2\pi\frac{1}{4}\right)^{2}=\frac{\pi}{8}$ in terms of our
proposed Ginzburg-Landau like energy.

The total singular energy is then given by

$\sum\limits_{i=1}^{m}\frac{(\pi-\alpha_{i}-k_{i}2\pi)^{2}}{2\alpha_{i}}+\frac{1}{2\pi}\left(\sum\limits_{i=1}^{n}\left(2\pi|l_{i}|\right)^{2}\right).$

Note that if we have $m$ boundary singularities of index $k_{i}$ and $n$
internal singularities of degrees $l_{j}$, the Poincaré-Hopf theorem implies
that

$\sum_{i=1}^{m}k_{i}+\sum_{j=1}^{n}l_{j}=\chi$

where $\chi$ is the Euler-Poincaré characteristic of the computational domain.

#### Optimality conditions

We first consider the case of balanced boundary singularities:
$\sum_{i=1}^{m}k_{i}=\chi$. It is the case for example if we have only four
boundary angles of $90^{\circ}$ that we associate with four 1/4 boundary
singularities with $\chi=1$.

######  Theorem 3.2.1 (Balanced boundary singularities).

If we have an optimal choice of boundary singularities in a case where the sum
of the index of the boundary equals the Euler-Poincaré characteristic of the
computational domain such as the Poincaré-Hopf theorem is fulfilled, the
limits angles $\alpha$ between which a boundary singularity of degree $k$ mst
be placed are given by the following relationship:

$\frac{1-\frac{1}{4}}{2}-\left(1+\frac{1}{8}\right)\frac{\alpha}{2\pi}\leq
k\leq\frac{1+\frac{1}{4}}{2}-\left(1-\frac{1}{8}\right)\frac{\alpha}{2\pi}.$
(15)

###### Proof.

See the appendix section 5.1.1. ∎

In particular, we can decide to place no singularity ($k_{f}=0$) when

$\frac{2\pi}{3}\leq\alpha_{f}\leq\frac{10\pi}{7}$

but the intervals overlap so we can decide to place a 1/4 singularity
($k_{f}=1/4$) when

$\frac{2\pi}{9}\leq\alpha_{f}\leq\frac{6\pi}{7}$

thus in the interval $2\pi/3\leq x\leq 6\pi/7$ we can place either a 1/4 or no
singularity. Note in particular that the heuristic classical limit angle of
$3\pi/4$ belongs to that interval.

For the two other cases, there is no overlaping between the intervals on which
place a singularity.

We now consider the case of rare boundary singularities:
$\sum_{i=1}^{m}k_{i}<\chi$. It is the case for example if we have three
boundary angles of $90^{\circ}$ that we associate with three 1/4 boundary
singularities with $\chi=1$.

######  Theorem 3.2.2 (Scarce boundary singularities).

If we have an optimal choice of boundary singularities in a case where the sum
of the index of the boundary singularities is lower than the Euler-Poincaré
characteristic of the computational domain, the limits angles $\alpha$ between
which a boundary singularity of degree $k$ must be placed are given by the
following relationship:

$\frac{1-\frac{1}{4}}{2}-\left(1-\frac{1}{8}\right)\frac{\alpha_{f}}{2\pi}\leq
k_{f}\leq\frac{1+\frac{1}{4}}{2}-\left(1-\frac{1}{8}\right)\frac{\alpha_{f}}{2\pi}.$
(16)

###### Proof.

See the appendix section 5.1.2. ∎

In particular, we can decide to place no singularity ($k_{f}=0$) when

$\frac{6\pi}{7}\leq\alpha_{f}\leq\frac{10\pi}{7}$

and we decide to place a 1/4 singularity ($k_{f}=1/4$) when

$\frac{2\pi}{7}\leq\alpha_{f}\leq\frac{6\pi}{7}.$

Note in particular that the heuristic classical limit angle of $3\pi/4$
belongs to that interval. This time the intervals follow each other without
overlapping:

$0\overset{k=\frac{2}{4}}{\longleftrightarrow}\frac{2\pi}{7}\overset{k=\frac{1}{4}}{\longleftrightarrow}\frac{6\pi}{7}\overset{k=0}{\longleftrightarrow}\frac{10\pi}{7}\overset{k=-\frac{1}{4}}{\longleftrightarrow}2\pi.$

We see that we place more easily positive singularity than negative ones.

We now consider the case of plentiful boundary singularities:
$\sum_{i=1}^{m}k_{i}>\chi$. It is the case for example if we have five
boundary angles of $90^{\circ}$ that we associate with five 1/4 boundary
singularities with $\chi=1$.

######  Theorem 3.2.3 (Excess of boundary singularities).

If we have an optimal choice of boundary singularities in a case where the sum
of the index of the boundary singularities is greater than the Euler-Poincaré
characteristic of the computational domain, the limits angles $\alpha$ between
which a boundary singularity of degree $k$ must be placed are given by the
following relationship:

$\frac{1-\frac{1}{4}}{2}-\left(1+\frac{1}{8}\right)\frac{\alpha_{f}}{2\pi}\leq
k_{f}\leq\frac{1+\frac{1}{4}}{2}-\left(1+\frac{1}{8}\right)\frac{\alpha_{f}}{2\pi}.$
(17)

###### Proof.

See the appendix section 5.1.3. ∎

In particular, we can decide to place no singularity ($k_{f}=0$) when

$\frac{2\pi}{3}\leq\alpha_{f}\leq\frac{10\pi}{9}$

and we decide to place a 1/4 singularity ($k_{f}=1/4$) when

$\frac{2\pi}{9}\leq\alpha_{f}\leq\frac{2\pi}{3}.$

Note in particular that the heuristic classical limit angle of $3\pi/4$
belongs to that interval. This time the intervals follow each other without
overlapping:

$0\overset{k=\frac{2}{4}}{\longleftrightarrow}\frac{2\pi}{9}\overset{k=\frac{1}{4}}{\longleftrightarrow}\frac{2\pi}{3}\overset{k=0}{\longleftrightarrow}\frac{10\pi}{9}\overset{k=-\frac{1}{4}}{\longleftrightarrow}\frac{14\pi}{9}\overset{k=-\frac{2}{4}}{\longleftrightarrow}2\pi.$

We see that we place more easily negative singularities than positive ones.

The conclusion of our calculations is that when we have enough boundary
singularities, the user can choose the type of the singularities on some
ranges; when we have not enough boundary singularities, the placement of
positive singularities occurs more often and when we have more than enough
boundary singularities, the placement of negative singularities occurs more
often. With the implementation, if we made a wrong supposition, choose the
another set of ranges will always be a better choice or an equivalent one if
we are in the situation where $\sum k_{i}=\chi$.

##  4 Results and Discussion

We will show and compare different results, as well as seeing that it is
sometimes necessary to add constraints to our basic formulation to obtain the
desired results. We will also give a new interpretation of the Poincaré-Hopf
theorem that is compatible with the optimizations proposed in the section 3.2.
This interpretation helps a user to properly apply our method without
necessarily knowing explicitly the configuration of the boundary
singularities.

###  4.1 Analysis of the numerical results

In figure 10, we can see that our formulation in $\mathbb{R}^{2}$ does not
result in a vector field in $S^{1}$ in the computed cross field and in the
corresponding vector field for which our algorithm has converged in 13 fast
iterations. Furthermore, we can easily move the singularities (see figure 11).

Fig 10: Cross field for a disk of radius 1 with four holes of radius $0.1$
with centers on a circle of radius $0.6$ at left and the corresponding vector
field at right.

Fig 11: Singularities placed via disks of radius equals to $0.01$ placed on
$(0,0)$, $(-0.2,0.2)$, $(0.2,-0.2)$ and $(-0.2,-0.2)$.

The method converges quickly and in never more than 30 iterations in the
tested cases. Furthermore, the Newton-Raphson scheme could be replaced with
methods better adapted to solve linear objective with quadratic constraints.
Even with the Newton-Raphson scheme the method can converge for any
topologically consistent configuration of the hole centers and sizes. On one
hand, larger holes will imply less rotations and thus less energy in the
computed cross fields but also less crosses in the computed field which will
thus imply quicker convergence as long as the configuration of singularity is
topologically coherent. On the other hand, large holes could imply more norm
disparity of the crosses on their outline and thus a worse estimation of the
degree on it. Furthermore, it is more difficult to extract a block-structured
quandrangular mesh from a domain perforated with large holes.

Among the results obtained, two are particularly worth mentioning.

First, for methods based on classical Ginzburg-Landau theories, the generation
of cross fields in a flat ellipse leads to singularities aligned on a single
line while the perforation method allows to keep a rectangular configuration
of singularities similar to that of the disc (see figure 12).

Fig 12: Ginzburg-Landau based singularity configuration (above) compared with
a perforated singularity configuration (below).

Second, we noticed that for a choice of singularities coherent and close to
the configuration prescribed by the Ginzburg-Landau functional (9), the
condition on a hole containing a singularity can be imposed on only one
triangle of the starting triangular mesh (see figure 13). If one departs too
far from the Ginzburg-Landau configuration, disparities of norms between the
crosses of the triangle may imply an imposition of unsatisfactory degree.

Fig 13: Generation of a cross field in which a singularity is placed on a
single triangle of the underlying mesh.

When we want to place a non-minimal number of singularities, we may have to
force the crosses on some holes to have a norm equal to 1 via additional
Lagrange multipliers (see the blue circle of the figure 14) to make the
solution converges towards the prescribed degrees. We exceptionally decided to
work with the domain $G$ (without perforated holes as singularities) rather
than $G_{\rho}$ in order to better visualize the different types of
singularities on figure 14. We also see that a large blue circle implies more
flexibility for the actual position of the placed singularity within the hole.
Note that we have here singularities on the corners boundary $\partial G$ and
that the Poincaré-Hopf theorem (3) is well respected thanks to those
singularities. Being very localized, the added constraints slow down the
method only moderately.

Fig 14: Forced singularities on a disk quarter. The constraint to obtain a
singularity -1/4 has been placed on the blue circle shown above. All crosses
are unitary on it. We see two compensatory 1/4 singularities appearing.

A flaw of the method is that the more we move from the solution associated
with the minimizers of the Ginzburg-Landau functional, the more we have
disparities on the norms of our crosses which can make the process of
extracting quadrangular meshes from it more difficult (see figure and 16). A
way to counteract this is to add to our method the non-linear penalty of (9)
via a linearization and a Newton-Raphson scheme that is performed until
convergence. All the computations implying this added terms have been
performed with $\epsilon=1/100$. There are two main drawbacks of this added
penalty. First, the modified method is far slower. Second, when we try to
place a holed singularity too far from the best location for the Ginzburg-
Landau functional, a zero norm singular cross of opposed degree sticks to it
and the natural zero norm singular cross of (9) appears (see figures 11 and
16).

Fig 15: Singularities obtained via the Ginzburg-Landau equation (9) on a disk
centered at the origin and of radius 1.

Fig 16: Singularities obtained via the Ginzburg-Landau equation (9) with some
large visible holes set to obtain singularities of degree 1/4 on a similar
disk as the one of figure 16.

Still there is at least one benefit with this added penalty term. With it, the
choices of singularities that do not respect the Poincaré-Hopf theorem are
permitted without divergence of the method. We can see on the figure 16 that a
wrong singularity has been placed on a circle centered on $(0.6,0.6)$ but two
opposite singularities appears on the hole to compensate it so as to fulfill
the Poincaré-Hopf theorem. We can also see that the wrongly placed
singularities on $(0.5,-0.5)$ have been canceled by sticking an opposite
singularity on it. Moreover, another singularity was placed on the inner
circle. Obviously, such cross fields are no longer regular.

To conclude this section, we propose some comments on the use of the method
and its possible extensions. On one hand, for a basic use one is invited to
set a size for the perforated holes similar to the mesh size of the underlying
triangular mesh and with one constraint on each boundary $\gamma_{i}$ of the
perforated holes. On the other hand, if you do not have a precise idea of the
configuration of the singularities, you can use a single global constraint
(14) as well as large circles as the boundaries of the singularities that you
may not perforate (as done in the figure 14) in order to more easily extract a
quadrangular cut. Note also that the choice of a Newton-Raphson scheme is
certainly not the most suited and that the alternating direction method of
multipliers (ADMM) may be more appropriate given our linear objective with
quadratic constraints. A very important application will be to compare energy
of different configuration (possibly non-minimal) such as to minimize them
(for example with positive and negative singularities). It is also important
to point out that our method does not allow us to impose boundary
singularities or to give information about them. To remedy this, we have
extended the Ginzburg-Landau theory in the section 3.2 and we will give a new
interpretation of the Poincaré-Hopf theorem in the section 4.2.

###  4.2 Generalized Poincaré-Hopf theorem

We now give a generalization of the Poincaré-Hopf theorem which is compatible
with the previously given optimizations of the boundary singularity degrees.
Given a geometry with boundary angles, this generalization allows to quickly
know what should be the sum of the degrees of the interior singularities.

The theoretical results obtained in the rich analysis of the energetic
Ginzurg-Landau functional do not apply to piecewise smooth domains (see [17]).
One way to get around this problem is to smooth the corners of a piecewise
smooth domain so as to be able to apply the Ginzburg-Landau theory.
Unfortunately, clear information about singularity degree for boundary
singularities which is seek in the literature (see [8]) can not be extracted
in this way. But in the section 3.2 we have created a new energy to get over
this fact.

On the other hand, multiple generalizations of the Poincaré-Hopf theorem have
been given in the literature (see [7], [8] and [3]). In [3] an extension of
the Poincaré-Hopf theorem to rational indices has first been developed.
However, these generalizations do not allow to have a priori information on
the configuration of singularities to be placed inside angular domains.
Indeed, the angles present on the boundary of such domains may or may not
correspond to singularities in the resulting mesh.

In this section, we propose a new formulation of the Poincaré-Hopf theorem for
N-directional fields to choose topologically coherent singularity
configuration to place inside angular domains.

In order to get information on internally placed singularity configurations,
we propose a strategy based on a valence term (illustrated on Figure 17) that
sums up the smallest oriented angles such that the object is mapped into an
equivalent object on each side of the discontinuity points of the boundary.

${\color[rgb]{0.25,0.46,0.02}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.46,0.02}\frac{3\pi}{4}}{\color[rgb]{0.25,0.46,0.02}\definecolor[named]{pgfstrokecolor}{rgb}{0.25,0.46,0.02}\leq\alpha<\pi}$${\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}\beta\
<\
}{\color[rgb]{0.82,0.01,0.11}\definecolor[named]{pgfstrokecolor}{rgb}{0.82,0.01,0.11}\frac{3\pi}{4}}$
Fig 17: On one hand, placing a $1/4$ singularity on a boundary angle $\alpha$
such as $3\pi/4\leq\alpha<\pi$ implies a lot of rotation (and therefore of
energy) in the computed cross field as seen on the section 3.2. On the other
hand, placing a 1/4 singularity on a small angle $\beta<3\pi/4$ implies much
less rotational energy in the cross field. The valence $V(X,\partial G)$ on
the two angles correspond to the angle between the successive oriented arrows
around the corners which is the minimal rotation that the crosses can perform
on those angles to keep their alignment to the boundary. In the case of
$\alpha$, this angle implies that no boundary singularity has been placed
since the arrow keeps pointing outside the domain. In the case of $\beta$ this
angle implies that a boundary singularity has been placed since the arrow
points outwards on one side of the angle and is parallel to the boundary on
the other side. Fig 18: In the first domain the purple orientation vector has
a sum of its smallest variations on the angles on the boundary such as to be
mapped on an equivalent cross equals to 0. Furthermore the sum of the indices
of the internal singularities equals 0 and since there are angles greater than
$\frac{3\pi}{2}$, there must be at least one singularity in the field
replacing them. For the second, the sum of the variations equals -1/2 and
there is 2 singularities -1/4 in the domain. For the third, the sum of the
variations equals 1/2 and there is 2 singularities 1/4 in the domain.

Let $\rho(M)$ be the set of the $M^{th}$ roots of the unity and $X$ be a
M-directional field defined on a piecewise smooth domain $G$ with a finite
number of angular points (jump points of the derivative of the
parametrization). We set

$\displaystyle\gamma:[0,1]\rightarrow\partial G,$
$\displaystyle\gamma(0)=\gamma(1)\text{ and }$ $\displaystyle\gamma\in
W^{1,\infty}([0,1]).$

Let $t_{j}$ be the angular points of $\partial G$ with $1\leq j\leq N-1$,
$t_{0}=0$ and $t_{N}=1$. Then we define the valence of a directional field $X$
on the boundary $\partial G$ as

$V(X,\partial
G)=\sum\limits_{j=1}^{N}\left(X^{+}(\gamma(t_{j}))-X^{-}(\gamma(t_{j}))\right)\left[\frac{2\pi}{\rho(M)}\right].$

Let $P_{1},P_{2},\ldots,P_{N}$ be the $N$ singularities of a field $X$ in the
interior of a domain $G$. Let $\phi$ be the angle between the tangent of the
domain boundary and an arbitrary fixed axis. For any piecewise smooth domains
we speculate that the following equality must hold :

$\sum\limits_{k=1}^{N}\text{index}(P_{k})=V(X,\partial
G)+\frac{1}{2\pi}\int_{\partial G\backslash\\{P_{k}\\}}d\phi.$ (18)

Therefore, the above figure of 18 presents a cross field with $V(X,\partial
G)=0$, the central figure a cross field with $V(X,\partial G)=-1/2$, and the
below figure a cross field with $V(X,\partial G)=1/2$. All the figures of this
article also corroborate the equation (18).

At this point, it is also important to note that the minimal number of
singularities on the first figure of 18 has been obtained even if the sum of
the internal singularities equals 0. With six boundary segments, we must place
singularities in order to end with 4-sided blocks. This configuration
minimizes the energy proposed in 2.6.1. These results corroborate the claims
made by Viertel and Osting [8] that Ginzburg-Landau based cross field
generation methods ensure the possibility to extract block-structured
quadrangular meshes from the computed cross fields. Furthermore, the presence
of inside singularities allow to “cancel” some of the angular singularities
that become regular boundary points from which originate three directions
towards the inside of the domain. It is cheaper to place one $1/4$ and one
$-1/4$ inside singularities than let the two angles close to $\pi$ at the left
of the domain be $1/4$ singularities and place two inside $-1/4$
singularities.

Let introduce an other cross field to corroborate equation (18) for domains
with Euler-Poincaré characteristics, $\chi$, different from 1.

Fig 19: Domain containing a centered regular hole and 6 singular holes (four
1/4 singular holes around the regular hole and two -1/4 between large angles).

Figure 19 is such that $V(X,\partial G)=1/2$ thanks to the external boundary
and $\int_{\partial G\backslash\\{x_{k}\\}}\frac{d\phi}{2\pi}=-1$ thanks to
the internal boundary. We must then have an internal singularity configuration
such that

$\sum\limits_{k=1}^{n}\text{index}(x_{k})=-1/2.$

However we have 4 flat angles that could be “canceled out” by 1/4
singularities. To complete those singularities four 1/4 singularities could be
placed to isolate the internal hole as in the configuration of holed
singularities proposed in 19.

In conclusion, the following algorithm could be use to place singularities
inside piecewise smooth closed domains. First, the equation (18) is computed.
Then, all the angles between $3\pi/4$ and $\pi$ not included are outlined and
the Euler-Poincaré characteristic is computed. Finally, a configurations of
holed singularities must be proposed. This configuration could be as to
“cancel out” the angles between $3\pi/4$ and $\pi$ not included. As our
methods does not allow to explicitly place boundary singularities, one must
verify that its configuration is such as to minimize the energy proposed in
3.2 for the boundary singularities while respecting the equation (18) and the
generalized Poincaré-Hopf theorem (3) with its imposed inside singularities.

###  4.3 Further thoughts

The minimizers of the Ginzburg-Landau energy

$E_{\epsilon}(u)=\frac{1}{2}\int_{G}|\nabla
u|^{2}+\frac{1}{4\epsilon^{2}}\int_{G}(|u|^{2}-1)^{2}$

are functions describing vector fields on surfaces.

The study of this energy has lead in its earlier developments to the
formulation of the problems used previously in this paper but the most notable
conclusion, as we have already seen, is that an equivalent common simpler
renormalized energy exists. This energy is about points and degrees of those
points on the surface of computation. We have seen that its computation
involves a scalar field $\Phi$. This scalar field $\Phi$ can be seen as the
harmonic conjugate of the phase $\theta$ of a complex valued vector field

$v(x,y)=\Phi(x,y)+i\theta(x,y).$

This is an interesting information for the computation of cross field since
the use of a function

$w(x,y)=H(x,y)+i\theta(x,y)$

is possible [31] to construct perpendicular fields $u$ and $v$ that form
together a cross field but the computation of $\Phi$ must be adapted to the
cross field generation. Notably, we must ensure alignment with the domain
boundaries. Let $H$ be such a field that we will give a way to calculate it
later. As it is the real part of a complex field, the scalar field $H$ is
directly linked with the size field of $w$. On the one hand, mathematicians
will be interested in the positions and degrees of a minimum number of
singularities that minimize the Ginzburg-Landau energy. On the other hand,
experts in quadrangular block-structured meshing will rather seek to impose
enough singularities in suited places and with suited degrees to produce a
cross field from which a quadrangular block-structured mesh respecting, for
example, a certain mesh size can be extracted. The computation of the field
$H$ is therefore very interesting in terms of size map of the subsequent
block-structured quadrangular mesh.

As far as cross-field generation is concerned, the scalar field $H$ is half of
a complex field from which a parameterization could be extracted (see [31]).
Moreover, we can place singularities in it at will while respecting the
Poincaré-Hopf theorem. $H(x,y)=\text{Re}(w(z))$ is solution of

$\left\\{\begin{array}[]{ccccc}\nabla^{2}H&=&2\pi\sum_{j=1}^{N}{k_{j}\over
4}\delta(x_{j},y_{j})&\text{in}&S\\\ \partial_{n}H&=&\kappa&\text{in}&\partial
S,\end{array}\right.$ (19)

where $\kappa$ is the local curvature of $\partial S$.

To construct the corresponding cross fields, ones can rely on a suited set of
cutting and the relation

$\nabla H=\nabla^{T}\theta$

but this remains a challenge.

A very interesting point of this new formulation is that it allows the
placement of boundary and interior singularities in one and the same way. This
allows a user to more easily and freely respect the Poincaré-Hopf theorem.

Here are an example of the computation of the field $H$ given the following
suitable set of singularities (the red points are $1/4$ singularities and the
green ones are $-1/4$ singularities such as $\sum_{j=1}^{N}k_{j}=4\chi(S)$).

Fig 20: Position of the singularities for the computation of the scalar field
$H$.

Fig 21: Computation of the scalar field $H$ from a linear Neumann problem
given the positions and orders of the deltas given in Figure 20 is given on
the left. On the right, we can see that the iso-values of the computed H are
perpendicular to the boundaries of the domain.

The problems that remain to be addressed are the extraction of a cross field,
and a size map based on a calculated $H$ field.

On this subject, note that size maps or local refinements already appear in
the literature which is a good thing for high-precision numerical simulation.

Jiang et al. [32] use them not as a goal but as a tool. They compute a
discrete metric on the input surface to obtain a cross field.

Lyon at al. [33] work directly on block-structured quadrangular meshes that
they locally refine. For this local refinement, they impose new singularities
i.e. new irregular vertices on the mesh. They are able to perform local
anisotropic refinement with few singularities. To do so they use a binary
program. They concentrate on local split of local edges or elements. Thus,
their method makes it possible to perform local surgeries on meshes already
generated where we recommend to take into account the size map as soon as the
cross-field is generated in order to directly produce the desired mesh.

## Conclusion

Most methods generating block-structured quadrangular meshes, with potential
user-based inputs, rely on the generation of smooth cross fields. In this
paper, we have broadened the theoretical context and given new methods to
generate such cross fields in the context of Ginzburg-Landau theory. In
particular, we have defined a new Ginzburg-Landau energy for the boundaries of
piecewise smooth domains containing singularities which was formally missing
in the scientific literature as noted by Viertel and Osting [8]. Moreover, we
have given a method to place singularities that does not correspond to those
of Ginzburg-Landau energy minimizers and offered a way to evaluate them in
terms of this energy. Practically, there are three ways to decrease the
Ginzburg-Landau energy while respecting the Poincaré-Hopf theorem: (i)
decrease the number of singularities, (ii) bring their degrees closer to zero,
or (iii) move them to make the cross field smoother.

What are the methods and theoretical elements presented in this article and
what can be said of their strengths and weaknesses. First, we have given a
method to generate a cross field with prescribed inner singularity
configurations. Said configurations consist of a set of holes with prescribed
degrees implying a different number of directions from the usual four that
meet on a cross. Numerically, the method corresponds to the optimization of a
linear objective function with very localized quadratic constraints. Another
constraint could be needed to ensure the degrees on these holes. Indeed, a
degree must be computed with crosses of a norm close to one. This method works
only with singularity configurations that respect the Poincaré-Hopf theorem.
Unfortunately, this method only allows to impose inner singularities and does
not give us clear information concerning the eventual boundary singularities.
Other generalizations of this method for singularity configurations that do
not respect the topology of the domain have been proposed but they are less
efficient, generate singular crosses in the field and unintended singularity
configurations. To help a user to respect the Poincaré-Hopf theorem we have
also proposed in this article new theoretical contributions.

In addition, we have given a way to evaluate the positions and degrees of such
holed singularities in terms of Ginzburg-Landau energy but without computing
an actual cross field which takes us away from the block-structured
quadrangular mesh generation. This evaluation gives the same Ginzburg-Landau
energy as the asymptotic minimizers of the Ginzburg-Landau functional if the
singularities configuration coincide and the radius of the holes tend towards
zero.

Hence, the proposed perforation method, as it both produces quality non-
singular cross fields and gives the option to choose the holed placements
mimicking inner singularities, could be a promising path to efficiently and
flexibly handle block-structured quadrangular mesh generation while keeping
quality criterion via Ginzburg-Landau energy.

We have ended the paper with a Neumann problem with zero-radius holes, similar
to the ones of the asymptotic Ginzburg-Landau energy for holed domains. In its
formulation both inner and boundary singularities are seen in the same way.
Furthermore, its input consists of configuration of both inner and boundary
singularities. This problem allows to construct a single-valued scalar field
strongly linked to the multivalued angles of a cross field. This field
corresponds to the harmonic conjugate of the phase field because when it is
added with i times the phase field, it gives a complex harmonic field almost
everywhere. Furthermore, the gradients of these two scalar fields are
perpendicular. Therefore, based on a given suited singularity configuration,
this scalar linear Neumann problem allows to construct a field of the harmonic
conjugates of the phases of a cross field. This problem may lead to a
particularly simple and informative way to compute cross field. Indeed, this
scalar field, as it correspond to the real part of a complex field, could be
strongly linked with the local size of the cross field.

Through the theoretical developments, we have constructed a Ginzburg-Landau
energy for the boundaries of piecewise smooth domains. Minimizing this energy
leads to new ranges of limit angles on which boundary singularities appear.
Besides, we have presented a new interpretation of the Poincaré-Hopf theorem.
This interpretation makes it possible to apprehend in a simpler way the total
degree of the configuration of the internal singularities. Although our main
method does not allow to place boundary singularities, those two contributions
help a user to gain crucial information about them.

Throughout this paper, we have used the energy of Ginzburg-Landau as a quality
criterion. We expect that the presented methods and theories can later be
adapted to other quality criteria, for instance, the adequacy of the future
block-structured quadrangular mesh to a prescribed size field.

We hope that our new contributions, unified under the Ginzburg-Landau theory,
could contribute to the construction of cross field with singularity
configuration focused on the needs of subsequent high performance numerical
simulations on block-structured quadrangular meshes.

## Acknowledgments

I would like to thank my dear friend Pierre de Wael for his generous help with
English writing. I would also like to thank Rémy Rodiac for his friendship and
his many clarifications.

##  5 Appendix

###  5.1 Demonstrations for the optimal ranges of boundary singularities

In the following demonstrations, we choose to add a $\pm 1/4$ singularity or
to change the order of a boundary singularity by adding $\pm 1/4$ to it. These
are concision choices but the argumentation holds true for other ones.

####  5.1.1 Demonstration of the theorem 3.2.1

###### Proof.

If we have an optimal choice of boundary singularities, we then should have by
optimality for every $i\in\\{1,\ldots,m\\}$,

$\sum_{i\neq
f}\left(\frac{\left(\frac{1}{2}-\frac{\alpha_{i}}{2\pi}-k_{i}\right)^{2}}{\frac{\alpha_{i}}{2\pi}}\right)+\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-\left(k_{f}\pm\frac{1}{4}\right)\right)^{2}}{\frac{\alpha_{f}}{2\pi}}+\sum_{j}|l_{j}|^{2}+\left(\frac{1}{4}\right)^{2}\\\
\geq\sum_{i}\left(\frac{\left(\frac{1}{2}-\frac{\alpha_{i}}{2\pi}-k_{i}\right)^{2}}{\frac{\alpha_{i}}{2\pi}}\right)+\sum_{j}|l_{j}|^{2},$

since if we have the optimal order of singularities, increase or decrease the
order of a boundary singularity (and add one compensating internal singularity
to respect the Poincaré-Hopf theorem) must imply more energy. Some
simplifications gives us the following inequality:

$\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-\left(k_{f}\pm\frac{1}{4}\right)\right)^{2}}{\frac{\alpha_{f}}{2\pi}}+\left(\frac{1}{4}\right)^{2}\geq\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}\right)^{2}}{\frac{\alpha_{f}}{2\pi}}.$

Developing this last inequality gives

$\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)^{2}-2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)\left(k_{f}\pm\frac{1}{4}\right)+\left(k_{f}\pm\frac{1}{4}\right)^{2}}{\frac{\alpha_{f}}{2\pi}}+\left(\frac{1}{4}\right)^{2}\geq\\\
\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)^{2}-2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)k_{f}+{k_{f}}^{2}}{\frac{\alpha_{f}}{2\pi}},$

which simplifies successively as

$\frac{\mp
2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)\frac{1}{4}+\left(k_{f}\pm\frac{1}{4}\right)^{2}}{\frac{\alpha_{f}}{2\pi}}+\left(\frac{1}{4}\right)^{2}\geq\frac{{k_{f}}^{2}}{\frac{\alpha_{f}}{2\pi}},$
$\frac{\mp 2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)\frac{1}{4}\pm
2k_{f}\frac{1}{4}+\left(\frac{1}{4}\right)^{2}}{\frac{\alpha_{f}}{2\pi}}+\left(\frac{1}{4}\right)^{2}\geq
0,$ $\frac{\left(\frac{1}{4}\right)^{2}\mp
2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}\right)\frac{1}{4}}{\frac{\alpha_{f}}{2\pi}}+\left(\frac{1}{4}\right)^{2}\geq
0,$

that we divide by $(1/4)^{2}$ and multiply by $\alpha_{f}/(2\pi)$ which gives

$1\mp
2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}\right)4+\frac{\alpha_{f}}{2\pi}\geq
0,$
$\frac{1}{8}\left(1+\frac{\alpha_{f}}{2\pi}\right)\geq\pm\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}\right),$
$-\frac{1}{8}\left(1+\frac{\alpha_{f}}{2\pi}\right)\leq\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}\right)\leq\frac{1}{8}\left(1+\frac{\alpha_{f}}{2\pi}\right),$
$-\frac{1}{2}+\frac{\alpha_{f}}{2\pi}-\frac{1}{8}\left(1+\frac{\alpha_{f}}{2\pi}\right)\leq-
k_{f}\leq-\frac{1}{2}+\frac{\alpha_{f}}{2\pi}+\frac{1}{8}\left(1+\frac{\alpha_{f}}{2\pi}\right),$
$\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-\frac{1}{8}\left(1+\frac{\alpha_{f}}{2\pi}\right)\leq
k_{f}\leq\frac{1}{2}-\frac{\alpha_{f}}{2\pi}+\frac{1}{8}\left(1+\frac{\alpha_{f}}{2\pi}\right),$
$\frac{1-\frac{1}{4}}{2}-\left(1+\frac{1}{8}\right)\frac{\alpha_{f}}{2\pi}\leq
k_{f}\leq\frac{1+\frac{1}{4}}{2}-\left(1-\frac{1}{8}\right)\frac{\alpha_{f}}{2\pi}.$

∎

####  5.1.2 Demonstration of the theorem 3.2.2

If we have an optimal choice of boundary singularities, we then should have by
optimality for every $i\in\\{1,\ldots,m\\}$,

$\sum_{i\neq
f}\left(\frac{\left(\frac{1}{2}-\frac{\alpha_{i}}{2\pi}-k_{i}\right)^{2}}{\frac{\alpha_{i}}{2\pi}}\right)+\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-\left(k_{f}\pm\frac{1}{4}\right)\right)^{2}}{\frac{\alpha_{f}}{2\pi}}+\sum_{j}|l_{j}|^{2}\\\
\leq\sum_{i}\left(\frac{\left(\frac{1}{2}-\frac{\alpha_{i}}{2\pi}-k_{i}\right)^{2}}{\frac{\alpha_{i}}{2\pi}}\right)+\sum_{j}|l_{j}|^{2}\pm\left(\frac{1}{4}\right)^{2},$

since we have boundary singularity degrees that sum up to less than $\chi$ but
we have choose the optimal choice of degrees by hypothesis, increase or
decrease the order of a boundary singularity must imply less energy than
adding or removing an internal one in order to closer respect the Poincaré-
Hopf theorem. Some simplifications gives us the following inequality:

$\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-\left(k_{f}\pm\frac{1}{4}\right)\right)^{2}}{\frac{\alpha_{f}}{2\pi}}\leq\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}\right)^{2}}{\frac{\alpha_{f}}{2\pi}}\pm\left(\frac{1}{4}\right)^{2}.$

Developing this last inequality gives

$\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)^{2}-2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)\left(k_{f}\pm\frac{1}{4}\right)+\left(k_{f}\pm\frac{1}{4}\right)^{2}}{\frac{\alpha_{f}}{2\pi}}\leq\\\
\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)^{2}-2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)k_{f}+{k_{f}}^{2}}{\frac{\alpha_{f}}{2\pi}}\pm\left(\frac{1}{4}\right)^{2},$

which simplifies successively as

$\frac{\mp
2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)\frac{1}{4}+\left(k_{f}\pm\frac{1}{4}\right)^{2}}{\frac{\alpha_{f}}{2\pi}}\leq\frac{{k_{f}}^{2}}{\frac{\alpha_{f}}{2\pi}}\pm\left(\frac{1}{4}\right)^{2},$
$\frac{\mp 2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)\frac{1}{4}\pm
2k_{f}\frac{1}{4}+\left(\frac{1}{4}\right)^{2}}{\frac{\alpha_{f}}{2\pi}}\leq\pm\left(\frac{1}{4}\right)^{2},$
$\frac{\left(\frac{1}{4}\right)^{2}\mp
2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}\right)\frac{1}{4}}{\frac{\alpha_{f}}{2\pi}}\leq\pm\left(\frac{1}{4}\right)^{2},$

that we divide by $(1/4)^{2}$ and multiply by $\alpha_{f}/(2\pi)$ which gives

$1\mp
2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}\right)4\leq\pm\frac{\alpha_{f}}{2\pi},$
$1\mp
2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}\right)4\mp\frac{\alpha_{f}}{2\pi}\leq
0,$ $\mp
2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}+\frac{1}{8}\frac{\alpha_{f}}{2\pi}\right)4+1\leq
0,$
$\frac{1}{8}\leq\pm\left(\frac{1}{2}+\frac{\alpha_{f}}{2\pi}\left(\frac{1}{8}-1\right)-k_{f}\right),$
$-\frac{1}{8}\leq\left(\frac{1}{2}+\frac{\alpha_{f}}{2\pi}\left(\frac{1}{8}-1\right)-k_{f}\right)\leq\frac{1}{8},$
$-\frac{1}{8}-\frac{1}{2}+\frac{\alpha_{f}}{2\pi}\left(1-\frac{1}{8}\right)\leq-
k_{f}\leq\frac{1}{8}-\frac{1}{2}+\frac{\alpha_{f}}{2\pi}\left(1-\frac{1}{8}\right),$
$\frac{1-\frac{1}{4}}{2}-\left(1-\frac{1}{8}\right)\frac{\alpha_{f}}{2\pi}\leq
k_{f}\leq\frac{1+\frac{1}{4}}{2}-\left(1-\frac{1}{8}\right)\frac{\alpha_{f}}{2\pi}.$

####  5.1.3 Demonstration of the theorem 3.2.3

If we have an optimal choice of boundary singularities, we then should have by
optimality for every $i\in\\{1,\ldots,m\\}$,

$\sum_{i\neq
f}\left(\frac{\left(\frac{1}{2}-\frac{\alpha_{i}}{2\pi}-k_{i}\right)^{2}}{\frac{\alpha_{i}}{2\pi}}\right)+\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-\left(k_{f}\mp\frac{1}{4}\right)\right)^{2}}{\frac{\alpha_{f}}{2\pi}}+\sum_{j}|l_{j}|^{2}\\\
\geq\sum_{i}\left(\frac{\left(\frac{1}{2}-\frac{\alpha_{i}}{2\pi}-k_{i}\right)^{2}}{\frac{\alpha_{i}}{2\pi}}\right)+\sum_{j}|l_{j}|^{2}\pm\left(\frac{1}{4}\right)^{2},$

since we have boundary singularity degrees that sum up to more that $\chi$ but
we have choose the optimal choice of degrees by hypothesis, increase or
decrease the order of a boundary singularity must imply more energy than
adding or removing an internal one in order to closer respect the Poincaré-
Hopf theorem. Some simplifications gives us the following inequality:

$\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-\left(k_{f}\mp\frac{1}{4}\right)\right)^{2}}{\frac{\alpha_{f}}{2\pi}}\geq\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}\right)^{2}}{\frac{\alpha_{f}}{2\pi}}\pm\left(\frac{1}{4}\right)^{2}.$

Developing this last inequality gives

$\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)^{2}-2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)\left(k_{f}\mp\frac{1}{4}\right)+\left(k_{f}\mp\frac{1}{4}\right)^{2}}{\frac{\alpha_{f}}{2\pi}}\geq\\\
\frac{\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)^{2}-2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)k_{f}+{k_{f}}^{2}}{\frac{\alpha_{f}}{2\pi}}\pm\left(\frac{1}{4}\right)^{2},$

which simplifies successively as

$\frac{\pm
2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)\frac{1}{4}+\left(k_{f}\pm\frac{1}{4}\right)^{2}}{\frac{\alpha_{f}}{2\pi}}\geq\frac{{k_{f}}^{2}}{\frac{\alpha_{f}}{2\pi}}\pm\left(\frac{1}{4}\right)^{2},$
$\frac{\pm 2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}\right)\frac{1}{4}\pm
2k_{f}\frac{1}{4}+\left(\frac{1}{4}\right)^{2}}{\frac{\alpha_{f}}{2\pi}}\geq\pm\left(\frac{1}{4}\right)^{2},$
$\frac{\left(\frac{1}{4}\right)^{2}\pm
2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}\right)\frac{1}{4}}{\frac{\alpha_{f}}{2\pi}}\geq\pm\left(\frac{1}{4}\right)^{2},$

that we divide by $(1/4)^{2}$ and multiply by $\alpha_{f}/(2\pi)$ which gives

$1\pm
2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}\right)4\geq\pm\frac{\alpha_{f}}{2\pi},$
$1\pm
2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}\right)4\mp\frac{\alpha_{f}}{2\pi}\geq
0,$ $\pm
2\left(\frac{1}{2}-\frac{\alpha_{f}}{2\pi}-k_{f}-\frac{1}{8}\frac{\alpha_{f}}{2\pi}\right)4+1\geq
0,$
$\frac{1}{8}\geq\pm\left(\frac{1}{2}+\frac{\alpha_{f}}{2\pi}\left(-\frac{1}{8}-1\right)-k_{f}\right),$
$-\frac{1}{8}\leq\left(\frac{1}{2}+\frac{\alpha_{f}}{2\pi}\left(-\frac{1}{8}-1\right)-k_{f}\right)\leq\frac{1}{8},$
$-\frac{1}{8}-\frac{1}{2}+\frac{\alpha_{f}}{2\pi}\left(1+\frac{1}{8}\right)\leq-
k_{f}\leq\frac{1}{8}-\frac{1}{2}+\frac{\alpha_{f}}{2\pi}\left(1+\frac{1}{8}\right),$
$\frac{1-\frac{1}{4}}{2}-\left(1+\frac{1}{8}\right)\frac{\alpha_{f}}{2\pi}\leq
k_{f}\leq\frac{1+\frac{1}{4}}{2}-\left(1+\frac{1}{8}\right)\frac{\alpha_{f}}{2\pi}.$

## References

  * 1. Ray N, Li WC, Lévy B, Sheffer A, Alliez P.  Periodic global parameterization.  ACM Transactions on Graphics. 2006;doi:10.1145/1183287.1183297.
  * 2. Palacios J, Zhang E.  Rotational Symmetry Field Design on Surfaces.  ACM Trans Graph. 2007;26(3):55–es.  doi:10.1145/1276377.1276446.
  * 3. Ray N, Vallet B, Li WC, Lévy B.  N-Symmetry Direction Field Design.  ACM Trans Graph. 2008;27(2).  doi:10.1145/1356682.1356683.
  * 4. Kälberer F, Nieser M, Polthier K; Wiley Online Library.  Quadcover-surface parameterization using branched coverings.  Computer graphics forum. 2007;26(3):375–384.
  * 5. Bommes D, Campen M, Ebke HC, Alliez P, Kobbelt L.  Integer-grid maps for reliable quad meshing.  ACM Transactions on Graphics (TOG). 2013;32(4):1–12.
  * 6. Campen M, Bommes D, Kobbelt L.  Quantized global parametrization.  Acm Transactions On Graphics (tog). 2015;34(6):1–12.
  * 7. Beaufort PA, Lambrechts J, Henrotte F, Geuzaine C, Remacle JF.  Computing cross fields - A PDE approach based on the Ginzburg-Landau theory.  Procedia Engineering. 2017;203:219–231.  doi:10.1016/j.proeng.2017.09.799.
  * 8. Viertel R, Osting B.  An Approach to Quad Meshing Based on Harmonic Cross-Valued Maps and the Ginzburg–Landau Theory.  SIAM Journal on Scientific Computing. 2019;41(1):A452–A479.  doi:10.1137/17M1142703.
  * 9. Vaxman A, Campen M, Diamanti O, Panozzo D, Bommes D, Hildebrandt K, et al.  Directional field synthesis, design, and processing.  In: Computer Graphics Forum; 2016.
  * 10. Bommes D, Lévy B, Pietroni N, Puppo E, Silva C, Tarini M, et al.  Quad Meshing.  In: Eurographics 2012 - State of the Art Reports; 2012.
  * 11. Bunin G.  A continuum theory for unstructured mesh generation in two dimensions.  Computer Aided Geometric Design. 2008;doi:10.1016/j.cagd.2007.05.002.
  * 12. Myles A, Zorin D.  Controlled-Distortion Constrained Global Parametrization.  ACM Trans Graph. 2013;32(4).  doi:10.1145/2461912.2461970.
  * 13. Ben-Chen M, Gotsman C, Bunin G.  Conformal flattening by curvature prescription and metric scaling.  In: Computer Graphics Forum; 2008.
  * 14. Sawhney R, Crane K.  Boundary First Flattening.  ACM Trans Graph. 2017;37(1):5:1–5:14.  doi:10.1145/3132705.
  * 15. Crane K, Desbrun M, Schröder P.  Trivial Connections on Discrete Surfaces.  Computer Graphics Forum (SGP). 2010;29(5):1525–1533.
  * 16. Hatcher A.  Algebraic topology.  Cambridge University Press; 2002.
  * 17. Bethuel F, Brezis H, Hélein F.  Ginzburg-Landau vortices.  Progress in Nonlinear Differential Equations and their Apllications. 1994;13.
  * 18. Macq A. Ginzburg-Landau et champs de croix : apports liés aux applications de maillage; 2018.  Available from: http://hdl.handle.net/2078.1/thesis:14770.
  * 19. Ginzburg VL, Landau LD.  On Superconductivity and Superfluidity: A Scientific Autobiography.  Springer; 2009.
  * 20. Kléman M.  Points, lignes, parois dans les fluides anisotropes et les solides cristallins.  Éditions de physique; 1977.
  * 21. Rubinstein J.  On the equilibrium position of Ginzburg Landau vortices.  ZAMP Zeitschrift für angewandte Mathematik und Physik. 1995;46(5):739–751.  doi:10.1007/BF00949077.
  * 22. Lassoued L, Mironescu P.  Ginzburg-Landau type energy with discontinuous constraint.  Journal d’Analyse Mathématique. 1999;77:1–26.  doi:10.1007/BF02791255.
  * 23. André N, Shafrir I.  Asymptotic behavior of minimizers for the Ginzburg-Landau functional with weight. Part II.  Archive for Rational Mechanics and Analysis. 1998;142(1):75–98.  doi:10.1007/s002050050084.
  * 24. Campen M, Zorin D.  Similarity maps and field-guided T-splines: a perfect couple.  ACM Transactions on Graphics. 2017;36:1–16.  doi:10.1145/3072959.3073647.
  * 25. Pietroni N, Puppo E, Marcias G, Roberto R, Cignoni P.  Tracing Field-Coherent Quad Layouts.  Computer Graphics Forum. 2016;35:485–496.  doi:10.1111/cgf.13045.
  * 26. Bommes D, Zimmer H, Kobbelt L.  Mixed-Integer Quadrangulation.  ACM Transactions on Graphics. 2009;28:article 77.  doi:10.1145/1576246.1531383.
  * 27. Pellenard B, Orbay G, Chen J, Sohan S, Kwok W, Tristano J.  QMCF: QMorph Cross Field-driven Quad-dominant Meshing Algorithm.  Procedia Engineering. 2014;82:338–350.  doi:10.1016/j.proeng.2014.10.395.
  * 28. Fang X, Bao H, Tong Y, Desbrun M, Huang J.  Quadrangulation through Morse-Parameterization Hybridization.  ACM Transactions on Graphics. 2018;37:1–15.  doi:10.1145/3197517.3201354.
  * 29. Crouzeix M, Raviart PA.  Conforming and nonconforming finite element methods for solving the stationary Stokes equations I.  ESAIM: Mathematical Modelling and Numerical Analysis. 1973;7(R3):33–75.
  * 30. Geuzaine C, Remacle JF.  Gmsh : a three-dimensional finite element mesh generator with built-in pre- and post- processing facilities.  International Journal for Numerical Methods in Engineering. 2009;79(11):1309–1331.
  * 31. Remacle JF, Henrotte F, Carrier Baudouin T, Geuzaine C, Béchet E, Mouton T, et al.  A frontal delaunay quad mesh generator using the L∞ norm.  In: Proceedings of the 20th International Meshing Roundtable, IMR 2011; 2011.
  * 32. Jiang T, Fang X, Huang J, Bao H, Tong Y, Desbrun M.  Frame field generation through metric customization.  In: ACM Transactions on Graphics; 2015.
  * 33. Lyon M, Bommes D, Kobbelt L.  Cost Minimizing Local Anisotropic Quad Mesh Refinement; 2020.