PDE-Based Block-Structured Meshing

Domain Partitioning via Cross Fields

Paper at a Glance

NoteReference

A PDE based approach to multi-domain partitioning and quadrilateral meshing — Kowalski, Ledoux, Frey (2012)

21st Int. Meshing Roundtable, Sandia Natl. Labs.

Download the original paper (PDF)

This page walks through every section of the paper in detail. The core idea: propagate boundary geometric information into the domain interior by solving a PDE, producing a smooth cross field whose topological structure (singularities and separatrices) partitions the domain into four-sided regions — each meshable by structured transfinite interpolation.

Try the interactive NACA airfoil demo

1. Introduction

“In numerous computational engineering applications, such as automobile crash simulations, structural mechanics, neutronics or fluid-structure interactions, quadrilateral meshes are preferred over triangular meshes.”

1.1 Why Quadrilateral Meshes?

The paper identifies four desirable properties of quad meshes:

  1. Minimize singularities — non-4-valent interior vertices (vertices that don’t have exactly 4 neighbors) should be rare
  2. Boundary alignment — several layers of quads aligned parallel to the boundary (boundary layers for CFD)
  3. Element quality — each quad should be as close to a square or rectangle as possible
  4. Size map conformance — elements should follow a prescribed size distribution

For multi-domain simulations (e.g., fluid-structure interaction), a fifth property is conformity along material interfaces.

1.2 Block-Structured Meshing

“Block-structured meshes have one main advantage over other kinds of meshes: they can theoretically adapt to any given domain while maintaining the good element quality associated with structured meshes.”

The key insight is that block-structured meshes combine the best of both worlds:

  • Structured meshes: high quality, efficient solvers (logical rectangular structure), hierarchical memory access — but can’t handle complex geometries
  • Unstructured meshes: handle any geometry — but lower quality, more singularities
  • Block-structured: unstructured arrangement of structured quad blocks. Each block is a four-sided region meshed by transfinite interpolation.

The challenge: how to automatically partition an arbitrary domain into four-sided regions? Previous attempts using the medial axis are numerically unstable and don’t produce four-sided regions (even a square gets triangular parts).

1.3 The Pipeline

The paper’s approach has four stages:

  1. Background triangulation \(T_h\) of domain \(\Omega\)
  2. Solve a PDE to propagate boundary orientation info into the interior, producing a smooth directionality (cross) field
  3. Analyze the topology of the cross field — find singularities and trace separatrices to partition the domain into four-sided regions
  4. Mesh each region with transfinite bilinear interpolation

2. Creation of a Directionality Field

“The main idea of our method consists in using a unit vector field prescribed on \(\partial\Omega\) that we propagate on the whole \(\Omega\).”

2.1 Definition of the Directionality of \(\Omega\)

A regular vertex \(P\) of a quad mesh is 4-valent. The four edges form two sets of opposite edges, each tracing a curve. The two curves intersect orthogonally at \(P\). Their tangent directions form a cross — the directionality of the mesh at \(P\).

A cross is a set of 4 unit vectors related by \(\pi/2\) rotations, fully described by a single angle \(\theta\):

\[ C_\theta = \left\{ \mathbf{v}_k = \left(\cos\!\left(\theta + \frac{k\pi}{2}\right), \sin\!\left(\theta + \frac{k\pi}{2}\right)\right)^T,\; 0 \le k \le 3 \right\}, \quad \theta \in \left[0, \frac{\pi}{2}\right). \]

The representation vector is \(\mathbf{u} = (\cos\theta, \sin\theta)^T\). This is not one of the cross arms — it encodes which of the infinitely many equivalent descriptions we pick.

Code 
viewof cross_theta_deg = Inputs.range([0, 89], {value: 25, step: 1, label: "θ (degrees)"})
Figure 1
Code 
{
  const w = 420, h = 420, cx = w/2, cy = h/2, r = 130;
  const theta = cross_theta_deg * Math.PI / 180;

  const svg = d3.create("svg").attr("viewBox", `0 0 ${w} ${h}`).style("width", "100%").style("max-width", w + "px");

  // Unit circle
  const circ = d3.range(0, 2*Math.PI + 0.01, 0.02).map(a => [cx + r*Math.cos(a), cy - r*Math.sin(a)]);
  svg.append("path").attr("d", d3.line()(circ)).attr("fill", "none").attr("stroke", "#ddd").attr("stroke-width", 1);

  // Axes
  svg.append("line").attr("x1", cx - r - 15).attr("y1", cy).attr("x2", cx + r + 15).attr("y2", cy).attr("stroke", "#eee");
  svg.append("line").attr("x1", cx).attr("y1", cy - r - 15).attr("x2", cx).attr("y2", cy + r + 15).attr("stroke", "#eee");

  // Cross arms
  const armColors = ["steelblue", "#e74c3c", "steelblue", "#e74c3c"];
  const armLabels = ["v0", "v1", "v2", "v3"];
  for (let k = 0; k < 4; k++) {
    const a = theta + k * Math.PI / 2;
    const dx = r * Math.cos(a), dy = -r * Math.sin(a);
    svg.append("line").attr("x1", cx).attr("y1", cy).attr("x2", cx + dx).attr("y2", cy + dy)
      .attr("stroke", armColors[k]).attr("stroke-width", 2.5).attr("stroke-linecap", "round");
    // Arrowhead
    const hl = 10, ha = 0.4;
    const ba = Math.atan2(-dy, dx) + Math.PI;
    svg.append("polygon")
      .attr("points", `${cx+dx},${cy+dy} ${cx+dx+hl*Math.cos(ba+ha)},${cy+dy-hl*Math.sin(ba+ha)} ${cx+dx+hl*Math.cos(ba-ha)},${cy+dy-hl*Math.sin(ba-ha)}`)
      .attr("fill", armColors[k]);
    svg.append("text").attr("x", cx + 1.15*dx).attr("y", cy + 1.15*dy)
      .attr("text-anchor", "middle").attr("dominant-baseline", "middle")
      .attr("font-size", 13).attr("fill", armColors[k]).text(armLabels[k]);
  }

  // Representation vector (orange, thicker)
  const repDx = r * 0.75 * Math.cos(theta), repDy = -r * 0.75 * Math.sin(theta);
  svg.append("line").attr("x1", cx).attr("y1", cy).attr("x2", cx + repDx).attr("y2", cy + repDy)
    .attr("stroke", "orange").attr("stroke-width", 4).attr("stroke-dasharray", "6,3");
  svg.append("text").attr("x", cx + repDx * 1.2).attr("y", cy + repDy * 1.2)
    .attr("text-anchor", "middle").attr("font-size", 14).attr("fill", "orange").attr("font-weight", "bold").text("u");

  // Angle arc
  if (theta > 0.02) {
    const arcPts = d3.range(0, theta + 0.01, 0.02).map(a => [cx + 35*Math.cos(a), cy - 35*Math.sin(a)]);
    svg.append("path").attr("d", d3.line()(arcPts)).attr("fill", "none").attr("stroke", "orange").attr("stroke-width", 2);
    svg.append("text").attr("x", cx + 48*Math.cos(theta/2)).attr("y", cy - 48*Math.sin(theta/2))
      .attr("text-anchor", "middle").attr("font-size", 14).attr("fill", "orange").text("θ");
  }

  svg.append("text").attr("x", cx).attr("y", 18).attr("text-anchor", "middle").attr("font-size", 15).attr("fill", "#333")
    .text(`Cross (θ = ${cross_theta_deg}°) with representation vector u`);

  return svg.node();
}

Key property: the cross field is independent of the reference axis. Rotating the coordinate system by \(\theta_0\) rotates each cross by \(-\theta_0\), and each representation vector by \(-4\theta_0\) (because of the 4-fold symmetry), leaving the problem globally unchanged.

“An important point to notice is that the cross field is not related to the reference axis used to compute angles.”

2.2 The PDE Formulation

The goal is to smoothly propagate boundary orientation information into the interior. This is analogous to the steady-state of a conduction (heat) problem:

\[ \frac{\partial \mathbf{u}}{\partial t} - \text{div}(k\nabla \mathbf{u}) = 0 \]

With \(k\) constant, at steady state we get the Laplace equation. Minimizing the Dirichlet energy gives us the smoothest possible field:

\[ \begin{cases} J(\mathbf{u}) = \int_\Omega |\nabla \mathbf{u}|^2\, dx, \\ \mathbf{u}(x) = \mathbf{u}_0(x) \quad \forall x \in \partial\Omega. \end{cases} \tag{Problem 3} \]

where \(\mathbf{u}\) is the representation vector field and \(\mathbf{u}_0\) is the Dirichlet boundary condition — the representation vector of the cross defined by tangent and normal directions at each boundary point.

Boundary conditions: At smooth boundary points, the cross is aligned with the tangent/normal. At \(C^0\) corners (sharp features), the representation vector is defined as the average of the two incident geometric edges’ representation vectors.

The Unit Norm Constraint

Problem 3 is a standard elliptic PDE — easy to solve. However, the solution may produce vectors with varying norms. Since only directions matter (not magnitudes), we add a non-linear constraint:

\[ \begin{cases} J(\mathbf{u}) = \int_\Omega |\nabla \mathbf{u}|^2\, dx, \\ \mathbf{u}(x) = \mathbf{u}_0(x) \quad \forall x \in \partial\Omega, \\ |\mathbf{u}(x)| = 1 \quad \forall x \in \Omega. \end{cases} \tag{Problem 4} \]

“Unlike Problem 2, Problem 4 is ill-posed, as a solution may not exist.”

The unit norm constraint prevents the field from having zeros. But the unconstrained solution may have zeros (where the field “wants” to vanish — these become singularities). The paper handles this through a two-step approach.

2.3 Solving the Linear Approximation

First, ignore the norm constraint and solve Problem 3. Splitting \(\mathbf{u} = \tilde{\mathbf{u}} + \mathbf{u}_0\), where \(\tilde{\mathbf{u}}\) vanishes on \(\partial\Omega\), the weak form is:

\[ \forall \mathbf{v} \in V, \quad \int_\Omega \nabla\tilde{\mathbf{u}} \nabla\mathbf{v}\, dx = -\int_\Omega \nabla\mathbf{u}_0 \nabla\mathbf{v}\, dx \]

with \(V = H^1_0(\Omega)\). The Lax-Milgram theorem guarantees existence and uniqueness.

Discretization: Galerkin FEM on the triangulation \(T_h\) with \(P_1\)-Lagrange elements (piecewise linear). This yields a symmetric positive-definite linear system \(Ax = b\), solved by conjugate gradient.

TipConnecting to the Demo

This is exactly what the “Solve linear” button does in the NACA airfoil demo — it solves the graph Laplacian (the discrete analog of the Dirichlet energy minimization) with Dirichlet conditions at constrained boundary vertices. The Gauss-Seidel iteration used there is equivalent to solving \(Ax = b\) iteratively.

2.4 Linearization of the Norm Constraint

The linear solution from Sec. 2.3 gives a good direction field, but vectors may not have unit norm. The paper refines this through an iterative process — alternating between a PDE solve and normalization.

At each step \(n\), given a current solution \(\mathbf{v}^n\) with \(|\mathbf{v}^n(p)| = 1\) at each mesh point \(p\), we add a linear constraint for step \(n+1\):

\[ \mathbf{v}^{n+1}(p) \cdot \mathbf{v}^n(p) = 1 \]

This is a linear approximation of \(|\mathbf{v}^{n+1}| = 1\) near the current solution. Using Lagrange multipliers \(\lambda_p\), this leads to a saddle-point system:

\[ \begin{pmatrix} A & C \\ C & 0 \end{pmatrix} \begin{pmatrix} u \\ \lambda \end{pmatrix} = \begin{pmatrix} b \\ 1 \end{pmatrix} \]

where \(C\) is a diagonal matrix encoding the constraint \(\mathbf{v}^{n+1} \cdot \mathbf{v}^n = 1\) at each interior point, with entries \(C(2p, 2p) = \alpha\) and \(C(2p+1, 2p+1) = \beta\) for \(\mathbf{v}^n(p) = (\alpha, \beta)\).

Algorithm: Solve the constrained system, then normalize all vectors to unit length. This guarantees \(|\mathbf{v}^{n+1}| \geq 1\) (the constraint ensures the component along \(\mathbf{v}^n\) is 1, so the full vector is at least unit length). Normalization then reduces \(J\). Repeat until convergence.

“The triangulation \(T_h\) that we use to solve the Problem 4 has no significative influence on the structure of the generated directionality field with regard to the convergence and the stability of the procedure.”

3. Domain Partitioning

“The general idea of our approach consists in analyzing the topology of this cross field. This means identifying its singularities, or zeroes, and connecting them, thus emphasizing its underlying structure.”

At this stage we have a smooth, unit-length cross field \(CF\) over the triangulation. The next step is to extract its topological skeleton — singularities connected by separatrices — which partitions the domain.

3.1 Singularities of a Directionality Field

Definition: A point \(x \in \Omega\) where \(\mathbf{v}(x) = 0\) is a singularity of the representation vector field. In practice, these are detected as triangles where the piecewise-linear interpolation of \(\mathbf{v}\) passes through zero.

The Poincare index measures the “topological charge” of a singularity:

\[ i_x = \frac{1}{2\pi} \oint_\gamma d\phi, \quad \text{with } \phi = \arctan\frac{v_1}{v_2}. \]

This counts the number of full rotations the vector field makes as you traverse a closed curve \(\gamma\) around the singularity.

Proposition 1: For a piecewise-linear vector field defined at triangle vertices with all non-zero values, all singularities are of first order — their index is either \(+1\) or \(-1\).

Code 
viewof singType = Inputs.radio(["index +1 (3-valent)", "index -1 (5-valent)"], {value: "index +1 (3-valent)", label: "Singularity type"})
Figure 2
Code 
{
  const w = 400, h = 400, cx = w/2, cy = h/2;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${w} ${h}`).style("width", "100%").style("max-width", w + "px");

  const isPos = singType.startsWith("index +1");
  const valence = isPos ? 3 : 5;
  const indexVal = isPos ? +1 : -1;

  // Draw crosses around the singularity
  const nCrosses = 12;
  const ringR = 140;
  const crossR = 22;

  // Singularity marker
  svg.append("circle").attr("cx", cx).attr("cy", cy).attr("r", 6).attr("fill", "#e74c3c");
  svg.append("text").attr("x", cx).attr("y", cy - 15).attr("text-anchor", "middle")
    .attr("font-size", 13).attr("fill", "#e74c3c").attr("font-weight", "bold")
    .text(`index = ${indexVal > 0 ? "+" : ""}${indexVal}`);

  for (let i = 0; i < nCrosses; i++) {
    const phi = 2 * Math.PI * i / nCrosses;
    const px = cx + ringR * Math.cos(phi);
    const py = cy - ringR * Math.sin(phi);

    // Cross orientation: for index +1, the cross rotates by +2pi going around once
    // In rep vector space: angle = indexVal * phi (for the cross field, index of rep vec = 4*index of cross)
    // For a +1 singularity of the cross field, rep vector index = +1, so rep rotates by 2pi
    // Cross itself rotates by pi/2 going around once
    const crossAngle = indexVal * phi / 4 + Math.PI / 8;

    for (let k = 0; k < 4; k++) {
      const a = crossAngle + k * Math.PI / 2;
      svg.append("line")
        .attr("x1", px + crossR * Math.cos(a)).attr("y1", py - crossR * Math.sin(a))
        .attr("x2", px - crossR * Math.cos(a)).attr("y2", py + crossR * Math.sin(a))
        .attr("stroke", k % 2 === 0 ? "steelblue" : "#888")
        .attr("stroke-width", 1.8).attr("stroke-linecap", "round");
    }
  }

  // Rotation arrow showing accumulated angle
  const arrowR = ringR + 35;
  const arcPts = d3.range(0, 1.8 * Math.PI, 0.02).map(a => [cx + arrowR * Math.cos(a), cy - arrowR * Math.sin(a)]);
  svg.append("path").attr("d", d3.line()(arcPts)).attr("fill", "none").attr("stroke", "#999").attr("stroke-width", 1).attr("stroke-dasharray", "4,3");

  svg.append("text").attr("x", cx).attr("y", h - 10).attr("text-anchor", "middle").attr("font-size", 14).attr("fill", "#333")
    .text(isPos ? "3 separatrices (valence 3)" : "5 separatrices (valence 5)");

  return svg.node();
}

The connection between representation vector field singularities and cross field structure:

  • Index \(+1\) of the vector field \(\Rightarrow\) the cross rotates by \(+\pi/2\) around the singularity \(\Rightarrow\) 3 separatrices emanate \(\Rightarrow\) 3-valent singularity in the quad mesh
  • Index \(-1\) of the vector field \(\Rightarrow\) the cross rotates by \(-\pi/2\) around the singularity \(\Rightarrow\) 5 separatrices emanate \(\Rightarrow\) 5-valent singularity in the quad mesh

3.2 Streamlines and Separatrices

A streamline of a vector field \(\mathbf{v}\) is a curve \(\gamma(s)\) whose tangent everywhere matches the field:

\[ \forall\, s \in [0,1],\quad \partial\gamma(s) \times \mathbf{v}(\gamma(s)) = 0 \]

For cross fields, a streamline can follow any of the 4 arms of the cross at each point:

\[ \forall\, s \in [0,1],\; \exists\, i \in [1..4],\quad \partial\gamma(s) \times \mathbf{v}_i(\gamma(s)) = 0 \]

A separatrix is a special streamline that ends at a singularity: \(\exists\, t \in [0,1],\; \gamma(t) = S_0\) where \(S_0\) is a singularity.

Numerical Computation of Separatrices

To trace a separatrix from singularity \(S_0\) located in triangle \(T_0\):

  1. Find exit points: For each edge of \(T_0\), check if a separatrix crosses it. At edge point \(P\), interpolate the cross field \(CF(P)\) from vertex values. A separatrix passes through \(P\) if the direction \(\overrightarrow{S_0 P}\) aligns with one of the 4 cross arms at \(P\):

\[ \exists\, k \in [0..3],\quad \frac{\mathbf{v}_k}{|\mathbf{v}_k|} = \frac{\mathbf{u}(P)}{|\mathbf{u}(P)|} \]

  1. Integrate through triangles: Once we have an exit point \(X_i\) on the boundary of triangle \(T\) with direction \(\mathbf{d}_i\), use Heun’s method (improved Euler) to find where the streamline exits \(T\):
    • The cross field is interpolated linearly from vertex representation vectors within \(T\)
    • The cross arm closest to the incoming direction \(\mathbf{d}_i\) is selected for consistency
    • The ODE \(Y = f(X)\) is integrated to find the next exit point \(X_{i+1}\) and direction \(\mathbf{d}_{i+1}\)
  2. Repeat until the separatrix hits a boundary (\(\partial\Omega\)) or another singularity.
TipConnecting to the Demo

The Streamlines visualization in the NACA airfoil demo traces cross field streamlines using exactly this approach — Euler integration through the mesh, selecting the closest cross arm at each step for direction continuity.

3.3 Definition of the Domain Partitioning

Combining the separatrices from all singularities (field singularities + geometric singularities at boundary corners) produces a partition of \(\Omega\).

“Every separatrix ends either on a singularity or on \(\partial\Omega\).”

Geometric singularities: The \(C^0\) corners of \(\partial\Omega\) are treated as additional singularities. Separatrices are traced from these corners into the domain using the same integration process.

Proposition 2: The set of separatrices of the geometric and field singularities leads to a partition of \(\Omega\) into 3- or 4-sided regions. The 3-sided regions necessarily contain a singularity.

Proof sketch: Every region is singularity-free (by construction, singularities are only at corners). Inside each region, the cross field is regular, so a parametrization \(f(u,v)\) exists. The boundaries inherit the cross field’s two orthogonal families of streamlines. A region can only be 3-sided if one edge degenerates to a point (a singularity where one family of streamlines converges).

Once we have four-sided regions, each is meshed independently via transfinite bilinear interpolation — a standard technique that maps the unit square \([0,1]^2\) to the curved quadrilateral region, producing a structured quad mesh within each block.

3.4 Deep Dive: Computing Streamline Directions

The pipeline demo in Section 6 traces streamlines through the cross field, but how exactly does the algorithm decide which direction to follow at each point? This section unpacks the math and the numerical scheme step by step.

Recall the streamline condition for a vector field (Definition 3 in the paper):

\[ \forall\, s \in [0,1], \quad \partial\gamma(s) \times \mathbf{v}(\gamma(s)) = 0 \]

The tangent of the curve \(\gamma\) at every point must be parallel to the field value at that point. For a cross field with four arms, we relax this to allow any of the four directions (Definition 4):

\[ \forall\, s \in [0,1],\; \exists\, i \in [1..4], \quad \partial\gamma(s) \times \mathbf{v}_i(\gamma(s)) = 0 \]

This means the streamline follows whichever cross arm is locally aligned with its current heading. The challenge is computing the field value at an arbitrary point inside a triangle — and maintaining directional continuity as the streamline crosses from one triangle to the next.

Step 1: Barycentric Interpolation of the Cross Field

Given a triangle with vertices \(S_1, S_2, S_3\) having representation vectors \(\mathbf{u}_1, \mathbf{u}_2, \mathbf{u}_3\), the field at any interior point \(P\) is interpolated using barycentric coordinates \((\lambda_1, \lambda_2, \lambda_3)\):

\[ \mathbf{u}(P) = \lambda_1 \mathbf{u}_1 + \lambda_2 \mathbf{u}_2 + \lambda_3 \mathbf{u}_3 \]

Since each representation vector encodes the cross angle \(\theta\) as \(\mathbf{u} = (\cos 4\theta, \sin 4\theta)\), the interpolation happens in representation space — not in angle space directly (averaging angles would fail near \(0/\pi\) wraparound). The cross angle at \(P\) is recovered as:

\[ \theta(P) = \frac{1}{4}\operatorname{atan2}\!\Big(\sum_k \lambda_k \sin 4\theta_k,\; \sum_k \lambda_k \cos 4\theta_k\Big) \]

This gives the four cross arms at \(P\) as:

\[ \mathbf{v}_k(P) = \left(\cos\!\left(\theta(P) + \frac{k\pi}{2}\right),\; \sin\!\left(\theta(P) + \frac{k\pi}{2}\right)\right), \quad k = 0,1,2,3 \]

Step 2: Arm Selection for Directional Continuity

At each integration step, we must pick one of the four arms to follow. The rule is simple: choose the arm whose direction is closest to the previous step’s direction \(\mathbf{d}_i\):

\[ k^* = \arg\max_{k \in \{0,1,2,3\}} \; \mathbf{v}_k(P) \cdot \mathbf{d}_i \]

This ensures the streamline doesn’t “jump” between arms as it crosses triangles — it smoothly follows one family of curves.

Try it below — click inside the triangle to place a probe point, and drag the green arrow to change the “previous direction.” Watch how the selected arm (highlighted) changes as you rotate the incoming direction.

Code 
viewof v0_angle_deg = Inputs.range([0, 89], {value: 10, step: 1, label: "Vertex A angle (°)"})
viewof v1_angle_deg = Inputs.range([0, 89], {value: 45, step: 1, label: "Vertex B angle (°)"})
viewof v2_angle_deg = Inputs.range([0, 89], {value: 75, step: 1, label: "Vertex C angle (°)"})
viewof prevDir_deg = Inputs.range([0, 359], {value: 30, step: 1, label: "Previous direction (°)"})
(a)
(b)
(c)
(d)
Figure 3
Code 
{
  const W = 560, H = 480;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`).style("width", "100%").style("max-width", W + "px");

  // Triangle vertices (large equilateral, centered)
  const triCx = W / 2, triCy = H / 2 + 20;
  const triR = 170;
  const verts = [
    { x: triCx + triR * Math.cos(Math.PI/2), y: triCy - triR * Math.sin(Math.PI/2), label: "A", angle: v0_angle_deg * Math.PI / 180 },
    { x: triCx + triR * Math.cos(Math.PI/2 + 2*Math.PI/3), y: triCy - triR * Math.sin(Math.PI/2 + 2*Math.PI/3), label: "B", angle: v1_angle_deg * Math.PI / 180 },
    { x: triCx + triR * Math.cos(Math.PI/2 + 4*Math.PI/3), y: triCy - triR * Math.sin(Math.PI/2 + 4*Math.PI/3), label: "C", angle: v2_angle_deg * Math.PI / 180 }
  ];

  // Draw triangle
  const triPath = `M${verts[0].x},${verts[0].y} L${verts[1].x},${verts[1].y} L${verts[2].x},${verts[2].y} Z`;
  svg.append("path").attr("d", triPath).attr("fill", "rgba(200,220,240,0.2)").attr("stroke", "#2c3e50").attr("stroke-width", 1.5);

  // Draw crosses at vertices
  const vertColors = ["steelblue", "#e74c3c", "#2c3e50"];
  const armLen = 28;
  verts.forEach((v, idx) => {
    for (let k = 0; k < 4; k++) {
      const a = v.angle + k * Math.PI / 2;
      svg.append("line")
        .attr("x1", v.x + armLen * Math.cos(a)).attr("y1", v.y - armLen * Math.sin(a))
        .attr("x2", v.x - armLen * Math.cos(a)).attr("y2", v.y + armLen * Math.sin(a))
        .attr("stroke", vertColors[idx]).attr("stroke-width", 2).attr("stroke-linecap", "round");
    }
    svg.append("text").attr("x", v.x + (v.x - triCx) * 0.22).attr("y", v.y + (v.y - triCy) * 0.22)
      .attr("text-anchor", "middle").attr("dominant-baseline", "middle")
      .attr("font-size", 15).attr("font-weight", "bold").attr("fill", vertColors[idx])
      .text(`${v.label} (${Math.round(v.angle * 180 / Math.PI)}°)`);
  });

  // Probe point at centroid (can be extended to click-to-place)
  const px = (verts[0].x + verts[1].x + verts[2].x) / 3;
  const py = (verts[0].y + verts[1].y + verts[2].y) / 3;

  // Barycentric coordinates
  const d = (verts[1].y - verts[2].y)*(verts[0].x - verts[2].x) + (verts[2].x - verts[1].x)*(verts[0].y - verts[2].y);
  const lam0 = ((verts[1].y - verts[2].y)*(px - verts[2].x) + (verts[2].x - verts[1].x)*(py - verts[2].y)) / d;
  const lam1 = ((verts[2].y - verts[0].y)*(px - verts[2].x) + (verts[0].x - verts[2].x)*(py - verts[2].y)) / d;
  const lam2 = 1 - lam0 - lam1;

  // Interpolate in (cos4θ, sin4θ) space
  const rc = lam0 * Math.cos(4*verts[0].angle) + lam1 * Math.cos(4*verts[1].angle) + lam2 * Math.cos(4*verts[2].angle);
  const rs = lam0 * Math.sin(4*verts[0].angle) + lam1 * Math.sin(4*verts[1].angle) + lam2 * Math.sin(4*verts[2].angle);
  const interpAngle = Math.atan2(rs, rc) / 4;

  // Draw probe point
  svg.append("circle").attr("cx", px).attr("cy", py).attr("r", 5).attr("fill", "orange");

  // Barycentric coordinate labels
  svg.append("text").attr("x", px + 10).attr("y", py - 25).attr("font-size", 11).attr("fill", "#666")
    .text(`λ = (${lam0.toFixed(2)}, ${lam1.toFixed(2)}, ${lam2.toFixed(2)})`);
  svg.append("text").attr("x", px + 10).attr("y", py - 10).attr("font-size", 11).attr("fill", "orange").attr("font-weight", "bold")
    .text(`θ(P) = ${(interpAngle * 180 / Math.PI).toFixed(1)}°`);

  // Draw interpolated cross arms at probe point
  const prevDirRad = prevDir_deg * Math.PI / 180;
  const probeArmLen = 40;
  let bestK = 0, bestDot = -Infinity;
  for (let k = 0; k < 4; k++) {
    const a = interpAngle + k * Math.PI / 2;
    const dot = Math.cos(a) * Math.cos(prevDirRad) + (-Math.sin(a)) * (-Math.sin(prevDirRad));
    if (dot > bestDot) { bestDot = dot; bestK = k; }
  }

  for (let k = 0; k < 4; k++) {
    const a = interpAngle + k * Math.PI / 2;
    const isSelected = (k === bestK);
    svg.append("line")
      .attr("x1", px + probeArmLen * Math.cos(a)).attr("y1", py - probeArmLen * Math.sin(a))
      .attr("x2", px - probeArmLen * Math.cos(a)).attr("y2", py + probeArmLen * Math.sin(a))
      .attr("stroke", isSelected ? "#27ae60" : "rgba(0,0,0,0.25)")
      .attr("stroke-width", isSelected ? 3.5 : 1.5)
      .attr("stroke-linecap", "round");
    if (isSelected) {
      // Arrowhead on selected arm
      const tipX = px + probeArmLen * Math.cos(a), tipY = py - probeArmLen * Math.sin(a);
      const ba = Math.atan2(probeArmLen * Math.sin(a), probeArmLen * Math.cos(a)) + Math.PI;
      const hl = 10, ha = 0.4;
      svg.append("polygon")
        .attr("points", `${tipX},${tipY} ${tipX+hl*Math.cos(ba+ha)},${tipY-hl*Math.sin(ba+ha)} ${tipX+hl*Math.cos(ba-ha)},${tipY-hl*Math.sin(ba-ha)}`)
        .attr("fill", "#27ae60");
      svg.append("text").attr("x", tipX + 12*Math.cos(a)).attr("y", tipY - 12*Math.sin(a))
        .attr("text-anchor", "middle").attr("font-size", 11).attr("fill", "#27ae60").attr("font-weight", "bold")
        .text("selected");
    }
  }

  // Draw previous direction arrow (green dashed)
  const pdLen = 55;
  const pdTipX = px + pdLen * Math.cos(prevDirRad), pdTipY = py - pdLen * Math.sin(prevDirRad);
  svg.append("line").attr("x1", px).attr("y1", py).attr("x2", pdTipX).attr("y2", pdTipY)
    .attr("stroke", "#27ae60").attr("stroke-width", 2.5).attr("stroke-dasharray", "6,4");
  const pdba = Math.atan2(pdLen*Math.sin(prevDirRad), pdLen*Math.cos(prevDirRad)) + Math.PI;
  svg.append("polygon")
    .attr("points", `${pdTipX},${pdTipY} ${pdTipX+8*Math.cos(pdba+0.4)},${pdTipY-8*Math.sin(pdba+0.4)} ${pdTipX+8*Math.cos(pdba-0.4)},${pdTipY-8*Math.sin(pdba-0.4)}`)
    .attr("fill", "#27ae60");
  svg.append("text").attr("x", pdTipX + 15*Math.cos(prevDirRad)).attr("y", pdTipY - 15*Math.sin(prevDirRad))
    .attr("text-anchor", "middle").attr("font-size", 11).attr("fill", "#27ae60")
    .text("d_prev");

  // Representation vector (orange dashed)
  const repLen = 35;
  svg.append("line").attr("x1", px).attr("y1", py)
    .attr("x2", px + repLen*Math.cos(interpAngle)).attr("y2", py - repLen*Math.sin(interpAngle))
    .attr("stroke", "orange").attr("stroke-width", 3).attr("stroke-dasharray", "5,3");

  // Title
  svg.append("text").attr("x", W/2).attr("y", 18).attr("text-anchor", "middle").attr("font-size", 14).attr("fill", "#333")
    .text("Cross field interpolation at probe point P (orange dot)");

  return svg.node();
}
TipWhy arm selection matters

Without arm selection, the streamline would randomly jump between the four cross directions at each step, producing a jagged, meaningless path. By always picking the arm closest to the previous heading, the algorithm commits to one “family” of streamlines and traces it smoothly. The orthogonal family is obtained by starting with a direction rotated by \(\pi/2\).

Step 3: Heun’s Method for Triangle-to-Triangle Integration

The paper uses Heun’s method (an improved Euler scheme) to trace streamlines through the triangulation. Given an entry point \(X_i\) on an edge of triangle \(T\) with direction \(\mathbf{d}_i\):

“We use the Heun’s method, an easy-to-use variation of a Runge-Kutta method, to integrate the ODE \(Y = f(X)\) over \(T\), starting from \(X_i\).”

Predictor (Euler step): Take a full step in the current direction to get a tentative exit point: \[ \tilde{X} = X_i + h \cdot \mathbf{v}_{\mathbf{d}_i}(X_i) \]

Corrector (average of endpoint slopes): Evaluate the field at both ends and average: \[ X_{i+1} = X_i + \frac{h}{2}\left[\mathbf{v}_{\mathbf{d}_i}(X_i) + \mathbf{v}_{\mathbf{d}_i}(\tilde{X})\right] \]

where \(h\) is chosen so the line segment exits the triangle. In practice, we find the intersection of the ray from \(X_i\) in direction \(\mathbf{d}_i\) with the opposite edge of \(T\).

The interactive below shows this process on a strip of connected triangles. Click Step to advance one triangle at a time, and toggle between Euler and Heun to see the improvement.

Click "Step" to trace through each triangle.

3.5 Deep Dive: Tracing Separatrices from Singularities

Separatrices are the structural skeleton of the cross field. They are special streamlines that start (or end) at singularities, and they form the boundaries of the block partition. Understanding how they are found and traced is essential.

“For each singularity \(S\), we build the separatrices, which are the streamlines of the cross field ending on \(S\).”

Finding Separatrix Exit Points (Eq. 11)

A singularity \(S_0\) sits inside a triangle \(T_0\) with vertices \(S_1, S_2, S_3\). The separatrices must exit \(T_0\) through its edges. The question is: where on each edge does a separatrix cross?

For any point \(P\) on an edge of \(T_0\), we can:

  1. Interpolate the cross field at \(P\) by linearly blending the vertex values (in representation space). This gives 4 unit vectors \(\mathbf{v}_0(P), \ldots, \mathbf{v}_3(P)\).

  2. Compute the direction from the singularity to \(P\): \(\mathbf{u}(P) = \overrightarrow{S_0 P} / |\overrightarrow{S_0 P}|\).

  3. Check for alignment: A separatrix passes through \(P\) if and only if one of the cross arms at \(P\) points directly toward (or away from) the singularity:

\[ \exists\, k \in [0..3], \quad \frac{\mathbf{v}_k(P)}{|\mathbf{v}_k(P)|} = \frac{\mathbf{u}(P)}{|\mathbf{u}(P)|} \]

As we sweep \(P\) along an edge, the cross field rotates smoothly while the direction \(\overrightarrow{S_0 P}\) also rotates. Each time an arm “catches up” to the S₀→P direction, we find an exit point. The number of catches depends on the singularity index:

  • Index +1: the cross makes an extra \(\pi/2\) rotation around \(S_0\) → 3 separatrices exit
  • Index -1: the cross is short by \(\pi/2\) → 5 separatrices exit
Code 
viewof sepIndexType = Inputs.radio(["index +1 (3 exits)", "index -1 (5 exits)"], {value: "index +1 (3 exits)", label: "Singularity type"})
viewof sepSampleDensity = Inputs.range([10, 50], {value: 25, step: 1, label: "Edge samples"})
(a)
(b)
Figure 4
Code 
{
  const W = 580, H = 520;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`).style("width", "100%").style("max-width", W + "px");

  const isPos = sepIndexType.startsWith("index +1");
  const expectedExits = isPos ? 3 : 5;

  // Triangle vertices (large)
  const cx = W/2, cy = H/2 + 10;
  const R = 200;
  const tri = [
    {x: cx, y: cy - R},
    {x: cx - R*Math.sin(2*Math.PI/3), y: cy + R*Math.cos(2*Math.PI/3)},
    {x: cx + R*Math.sin(2*Math.PI/3), y: cy + R*Math.cos(2*Math.PI/3)}
  ];

  // Vertex cross angles chosen to produce the desired singularity index.
  // For index +1: angles that sum to +π/2 of angle defect
  // For index -1: angles that sum to -π/2 of angle defect
  const vertAngles = isPos
    ? [0.0, Math.PI/6, Math.PI/3 + Math.PI/8]
    : [0.0, -Math.PI/8, Math.PI/6];

  // Singularity at centroid
  const sx = (tri[0].x + tri[1].x + tri[2].x) / 3;
  const sy = (tri[0].y + tri[1].y + tri[2].y) / 3;

  // Draw triangle
  svg.append("path")
    .attr("d", `M${tri[0].x},${tri[0].y} L${tri[1].x},${tri[1].y} L${tri[2].x},${tri[2].y} Z`)
    .attr("fill", "rgba(200,220,240,0.12)").attr("stroke", "#2c3e50").attr("stroke-width", 1.5);

  // Draw vertex crosses
  const armLen = 22;
  const vColors = ["steelblue", "#e74c3c", "#8e44ad"];
  for (let vi = 0; vi < 3; vi++) {
    for (let k = 0; k < 4; k++) {
      const a = vertAngles[vi] + k * Math.PI / 2;
      svg.append("line")
        .attr("x1", tri[vi].x + armLen*Math.cos(a)).attr("y1", tri[vi].y - armLen*Math.sin(a))
        .attr("x2", tri[vi].x - armLen*Math.cos(a)).attr("y2", tri[vi].y + armLen*Math.sin(a))
        .attr("stroke", vColors[vi]).attr("stroke-width", 1.8).attr("stroke-linecap", "round");
    }
  }

  // Singularity marker
  svg.append("circle").attr("cx", sx).attr("cy", sy).attr("r", 7).attr("fill", "#e74c3c");
  svg.append("text").attr("x", sx).attr("y", sy - 14).attr("text-anchor", "middle")
    .attr("font-size", 12).attr("fill", "#e74c3c").attr("font-weight", "bold").text("S₀");

  // Sample each edge and check for alignment
  const edges = [[0,1],[1,2],[2,0]];
  const exitPoints = [];
  const allSamples = [];
  const nSamp = sepSampleDensity;
  const alignThreshold = Math.PI / nSamp * 2.0;

  for (const [ea, eb] of edges) {
    for (let si = 1; si < nSamp; si++) {
      const t = si / nSamp;
      const px = tri[ea].x + t * (tri[eb].x - tri[ea].x);
      const py = tri[ea].y + t * (tri[eb].y - tri[ea].y);

      // Barycentric coords for this point
      const dd = (tri[1].y-tri[2].y)*(tri[0].x-tri[2].x)+(tri[2].x-tri[1].x)*(tri[0].y-tri[2].y);
      const l0 = ((tri[1].y-tri[2].y)*(px-tri[2].x)+(tri[2].x-tri[1].x)*(py-tri[2].y))/dd;
      const l1 = ((tri[2].y-tri[0].y)*(px-tri[2].x)+(tri[0].x-tri[2].x)*(py-tri[2].y))/dd;
      const l2 = 1-l0-l1;

      // Interpolate cross field
      const rc = l0*Math.cos(4*vertAngles[0]) + l1*Math.cos(4*vertAngles[1]) + l2*Math.cos(4*vertAngles[2]);
      const rs = l0*Math.sin(4*vertAngles[0]) + l1*Math.sin(4*vertAngles[1]) + l2*Math.sin(4*vertAngles[2]);
      const baseAngle = Math.atan2(rs, rc) / 4;

      // Direction S0 -> P
      const dx = px - sx, dy = py - sy;
      const dirToP = Math.atan2(-dy, dx); // SVG y is flipped

      // Check each arm for alignment
      let aligned = false;
      let alignedArm = -1;
      for (let k = 0; k < 4; k++) {
        const arm = baseAngle + k * Math.PI / 2;
        let diff = arm - dirToP;
        diff = diff - Math.round(diff / (2*Math.PI)) * 2*Math.PI;
        if (Math.abs(diff) < alignThreshold) {
          aligned = true;
          alignedArm = k;
          break;
        }
      }

      allSamples.push({px, py, baseAngle, dirToP, aligned, alignedArm});

      if (aligned) {
        // Check it's not a duplicate
        let dup = false;
        for (const ep of exitPoints) {
          if (Math.sqrt((px-ep.px)**2+(py-ep.py)**2) < 20) { dup = true; break; }
        }
        if (!dup) exitPoints.push({px, py, dir: baseAngle + alignedArm * Math.PI/2});
      }
    }
  }

  // Draw sample points and mini-crosses
  const miniArm = 10;
  for (const s of allSamples) {
    // Direction line from S0 to P (thin gray)
    svg.append("line").attr("x1", sx).attr("y1", sy).attr("x2", s.px).attr("y2", s.py)
      .attr("stroke", "rgba(0,0,0,0.04)").attr("stroke-width", 0.5);

    // Sample dot
    svg.append("circle").attr("cx", s.px).attr("cy", s.py).attr("r", s.aligned ? 4 : 2)
      .attr("fill", s.aligned ? "#27ae60" : "rgba(0,0,0,0.15)");

    // Mini cross at sample point (only every 3rd point to avoid clutter)
    if (allSamples.indexOf(s) % 3 === 0) {
      for (let k = 0; k < 4; k++) {
        const a = s.baseAngle + k * Math.PI / 2;
        const isAlignedArm = s.aligned && (k === s.alignedArm);
        svg.append("line")
          .attr("x1", s.px + miniArm*Math.cos(a)).attr("y1", s.py - miniArm*Math.sin(a))
          .attr("x2", s.px - miniArm*Math.cos(a)).attr("y2", s.py + miniArm*Math.sin(a))
          .attr("stroke", isAlignedArm ? "#27ae60" : "rgba(0,0,0,0.12)")
          .attr("stroke-width", isAlignedArm ? 2.5 : 0.8);
      }
    }
  }

  // Highlight exit points
  for (const ep of exitPoints) {
    svg.append("circle").attr("cx", ep.px).attr("cy", ep.py).attr("r", 7)
      .attr("fill", "none").attr("stroke", "#27ae60").attr("stroke-width", 2.5);
    // Direction arrow from exit point
    const arrowLen = 30;
    const tipX = ep.px + arrowLen * Math.cos(ep.dir);
    const tipY = ep.py - arrowLen * Math.sin(ep.dir);
    svg.append("line").attr("x1", ep.px).attr("y1", ep.py).attr("x2", tipX).attr("y2", tipY)
      .attr("stroke", "#27ae60").attr("stroke-width", 2.5).attr("stroke-linecap", "round");
    // Small arrowhead
    const ba = Math.atan2(arrowLen*Math.sin(ep.dir), arrowLen*Math.cos(ep.dir)) + Math.PI;
    svg.append("polygon")
      .attr("points", `${tipX},${tipY} ${tipX+7*Math.cos(ba+0.4)},${tipY-7*Math.sin(ba+0.4)} ${tipX+7*Math.cos(ba-0.4)},${tipY-7*Math.sin(ba-0.4)}`)
      .attr("fill", "#27ae60");
  }

  // Label
  svg.append("text").attr("x", W/2).attr("y", H - 8).attr("text-anchor", "middle")
    .attr("font-size", 13).attr("fill", "#333")
    .text(`Found ${exitPoints.length} exit points (expected ${expectedExits}) — green circles show separatrix exits`);

  svg.append("text").attr("x", W/2).attr("y", 18).attr("text-anchor", "middle")
    .attr("font-size", 14).attr("fill", "#333")
    .text(`Separatrix exit detection: ${sepIndexType}`);

  return svg.node();
}

Tracing Outward from Exit Points

Once the exit points are found, the algorithm from Section 3.4 takes over: starting from each exit point with its known direction, we integrate outward through successive triangles using Heun’s method. The separatrix continues until it hits:

  • The domain boundary \(\partial\Omega\)
  • Another singularity
  • Another separatrix (forming a T-junction)

Each traced separatrix becomes an edge in the block partition graph.

The interactive below shows separatrices growing outward from a singularity embedded in a small mesh patch. Each separatrix is colored distinctly.

Click "Trace" to grow separatrices from the singularity.

3.6 Deep Dive: From Separatrices to Block Partition

Once all separatrices are traced — from both field singularities (interior points where the representation vector vanishes) and geometric singularities (\(C^0\) corners of \(\partial\Omega\)) — they divide the domain into regions.

Geometric Singularities

“In order to get four-sided regions, we extend the set of representation vector field singularities by adding a set of geometric singularities, which are defined as the \(C^0\) corners on \(\partial\Omega\).”

A sharp corner on the boundary (like the corners of an L-shape, or the trailing edge of an airfoil) acts as a singularity source. The representation vector at a corner is defined as the average of the representation vectors of the two incident boundary edges. Separatrices from geometric corners are traced inward using the same Heun integration.

Proposition 2: The Partition Theorem

Proposition 2 (Kowalski et al.): The set of separatrices of the geometric and field singularities leads to a partition of \(\Omega\) into 3-sided or 4-sided regions. The 3-sided regions necessarily contain a singularity.

Why only 3 or 4 sides? Inside each region (away from singularities), the cross field is smooth and regular. This means there exist two orthogonal families of curves — the \(u\)-family and the \(v\)-family — that tile the region like a coordinate grid. The boundary of each region consists of alternating segments from these two families:

  • A 4-sided region has two sides from the \(u\)-family and two from the \(v\)-family — like a curved rectangle. This is the ideal case for transfinite interpolation.
  • A 3-sided region occurs when one “side” degenerates to a single point — a singularity where all streamlines of one family converge. The angle at this degenerate vertex is at most \(\pi/2\).
Code 
{
  const W = 600, H = 280;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`).style("width", "100%").style("max-width", W + "px");

  // 4-sided region (left)
  const leftCx = 150, topY = 40, botY = 240;
  // Draw curved quad
  const quadPts = [
    {x: 60, y: 60}, {x: 240, y: 50}, {x: 260, y: 230}, {x: 50, y: 240}
  ];
  svg.append("path")
    .attr("d", `M${quadPts[0].x},${quadPts[0].y} C${quadPts[0].x+60},${quadPts[0].y-10} ${quadPts[1].x-60},${quadPts[1].y-10} ${quadPts[1].x},${quadPts[1].y} L${quadPts[2].x},${quadPts[2].y} C${quadPts[2].x-60},${quadPts[2].y+5} ${quadPts[3].x+60},${quadPts[3].y+5} ${quadPts[3].x},${quadPts[3].y} Z`)
    .attr("fill", "rgba(46,204,113,0.15)").attr("stroke", "#2c3e50").attr("stroke-width", 2);

  // Draw u-family curves (horizontal-ish)
  for (let i = 1; i <= 3; i++) {
    const t = i/4;
    const lx = quadPts[0].x + t*(quadPts[3].x-quadPts[0].x);
    const ly = quadPts[0].y + t*(quadPts[3].y-quadPts[0].y);
    const rx = quadPts[1].x + t*(quadPts[2].x-quadPts[1].x);
    const ry = quadPts[1].y + t*(quadPts[2].y-quadPts[1].y);
    svg.append("line").attr("x1",lx).attr("y1",ly).attr("x2",rx).attr("y2",ry)
      .attr("stroke","steelblue").attr("stroke-width",1).attr("stroke-dasharray","4,3").attr("opacity",0.6);
  }
  // Draw v-family curves (vertical-ish)
  for (let i = 1; i <= 3; i++) {
    const t = i/4;
    const tx = quadPts[0].x + t*(quadPts[1].x-quadPts[0].x);
    const ty = quadPts[0].y + t*(quadPts[1].y-quadPts[0].y);
    const bx = quadPts[3].x + t*(quadPts[2].x-quadPts[3].x);
    const by = quadPts[3].y + t*(quadPts[2].y-quadPts[3].y);
    svg.append("line").attr("x1",tx).attr("y1",ty).attr("x2",bx).attr("y2",by)
      .attr("stroke","#e74c3c").attr("stroke-width",1).attr("stroke-dasharray","4,3").attr("opacity",0.6);
  }

  svg.append("text").attr("x",150).attr("y",270).attr("text-anchor","middle").attr("font-size",13).attr("fill","#333")
    .text("4-sided region (ideal)");

  // Corner dots
  for (const p of quadPts) {
    svg.append("circle").attr("cx",p.x).attr("cy",p.y).attr("r",4).attr("fill","#2c3e50");
  }

  // Labels
  svg.append("text").attr("x",150).attr("y",35).attr("text-anchor","middle").attr("font-size",11).attr("fill","steelblue").text("u-family");
  svg.append("text").attr("x",275).attr("y",145).attr("font-size",11).attr("fill","#e74c3c").text("v-family");

  // 3-sided region (right)
  const ox = 330;
  // Triangle with one degenerate vertex (top = singularity)
  const triPts = [
    {x: ox + 150, y: 50},   // singularity (degenerate vertex)
    {x: ox + 260, y: 235},
    {x: ox + 40, y: 230}
  ];

  svg.append("path")
    .attr("d", `M${triPts[0].x},${triPts[0].y} C${triPts[0].x+50},${triPts[0].y+60} ${triPts[1].x-20},${triPts[1].y-60} ${triPts[1].x},${triPts[1].y} L${triPts[2].x},${triPts[2].y} C${triPts[2].x+20},${triPts[2].y-60} ${triPts[0].x-50},${triPts[0].y+60} ${triPts[0].x},${triPts[0].y}`)
    .attr("fill", "rgba(231,76,60,0.1)").attr("stroke", "#2c3e50").attr("stroke-width", 2);

  // Streamlines converging to the singularity
  for (let i = 1; i <= 4; i++) {
    const t = i/5;
    const bx = triPts[2].x + t*(triPts[1].x - triPts[2].x);
    const by = triPts[2].y + t*(triPts[1].y - triPts[2].y);
    svg.append("line").attr("x1",triPts[0].x).attr("y1",triPts[0].y).attr("x2",bx).attr("y2",by)
      .attr("stroke","steelblue").attr("stroke-width",1).attr("stroke-dasharray","4,3").attr("opacity",0.5);
  }

  // Singularity marker
  svg.append("circle").attr("cx",triPts[0].x).attr("cy",triPts[0].y).attr("r",6).attr("fill","#e74c3c");
  svg.append("text").attr("x",triPts[0].x).attr("y",triPts[0].y-12).attr("text-anchor","middle")
    .attr("font-size",11).attr("fill","#e74c3c").attr("font-weight","bold").text("singularity");

  // Other corners
  svg.append("circle").attr("cx",triPts[1].x).attr("cy",triPts[1].y).attr("r",4).attr("fill","#2c3e50");
  svg.append("circle").attr("cx",triPts[2].x).attr("cy",triPts[2].y).attr("r",4).attr("fill","#2c3e50");

  svg.append("text").attr("x",ox+150).attr("y",270).attr("text-anchor","middle").attr("font-size",13).attr("fill","#333")
    .text("3-sided region (degenerate vertex)");

  return svg.node();
}
Figure 5

Interactive: Notched Domain Partition

This domain — a rectangle with two semicircular notches — is a classic test case from Kowalski et al. The curved concavities force the cross field to develop interior singularities whose separatrices partition the domain into quad-meshable blocks. Step through each stage to see the full pipeline.

Click "Mesh" to generate the notched domain triangulation.
NoteWhat you should notice
  • The curved concavities force field singularities to appear in the interior — unlike polygonal domains where corners drive the topology
  • Singularities appear near the narrowest part of the domain (the “waist”) where the cross field must rotate significantly to stay aligned with the curved boundaries
  • Each colored region is a candidate for structured meshing — 4-sided regions can be directly meshed by transfinite interpolation

3.7 Deep Dive: Block Merging and Transfinite Interpolation

The separatrix-based partition from Section 3.6 may not immediately produce all 4-sided regions. Some regions may be 3-sided (with a degenerate vertex at a singularity), and some partitions may have more blocks than necessary. Block merging is a post-processing step that simplifies the partition.

Why Merge Blocks?

  1. 3-sided regions have a degenerate edge (one “side” is a single point — the singularity). Standard transfinite interpolation requires 4 distinct boundary curves.
  2. Excessive blocks create unnecessary complexity: each block interface must be discretized consistently, and more blocks mean more interface constraints.

The Merging Algorithm

When a 3-sided region \(R_3\) shares a separatrix with an adjacent 4-sided region \(R_4\), and they meet at a valence-3 singularity, we can remove the shared separatrix to merge them into a new 4-sided region:

while any 3-sided region R3 exists:
    find adjacent 4-sided region R4 sharing separatrix S at singularity V
    remove S from the partition
    merge R3 and R4 into new 4-sided region R_new
    update partition graph

The topological argument: removing a separatrix between a 3-sided and 4-sided region that share a vertex eliminates the degenerate edge and produces a valid 4-sided region. The combined region inherits 3 sides from \(R_4\) and 2 from \(R_3\) — but the shared separatrix accounts for one side from each, so the result has \(4 + 3 - 2 = 4\) sides (accounting for the singularity vertex no longer being a degenerate edge but a regular boundary vertex).

Note

This merging step is not explicitly described in Kowalski et al. 2012, but it is standard practice in block-structured meshing pipelines. The paper’s Proposition 2 guarantees that the initial partition produces only 3- or 4-sided regions, and the merging step eliminates all 3-sided ones.

Code 
viewof mergeStep = Inputs.range([0, 3], {value: 0, step: 1, label: "Merge step"})
Figure 6
Code 
{
  const W = 560, H = 400;
  const svg = d3.create("svg").attr("viewBox", `0 0 ${W} ${H}`).style("width", "100%").style("max-width", W + "px");

  // Pre-defined block partition with 3- and 4-sided regions
  // Simple example: a rectangular domain split by separatrices from 2 singularities
  const regions = [
    { // Region 0: 4-sided (top-left)
      pts: [{x:30,y:30},{x:270,y:30},{x:250,y:160},{x:30,y:170}],
      sides: 4, color: 'rgba(52,152,219,0.25)', label: '4'
    },
    { // Region 1: 3-sided (top-center triangle)
      pts: [{x:270,y:30},{x:400,y:30},{x:250,y:160}],
      sides: 3, color: 'rgba(231,76,60,0.2)', label: '3'
    },
    { // Region 2: 4-sided (right)
      pts: [{x:400,y:30},{x:530,y:30},{x:530,y:190},{x:320,y:180}],
      sides: 4, color: 'rgba(46,204,113,0.25)', label: '4'
    },
    { // Region 3: 3-sided (center triangle)
      pts: [{x:250,y:160},{x:320,y:180},{x:280,y:260}],
      sides: 3, color: 'rgba(155,89,182,0.2)', label: '3'
    },
    { // Region 4: 4-sided (bottom-left)
      pts: [{x:30,y:170},{x:250,y:160},{x:280,y:260},{x:30,y:370}],
      sides: 4, color: 'rgba(241,196,15,0.25)', label: '4'
    },
    { // Region 5: 4-sided (bottom-right)
      pts: [{x:280,y:260},{x:320,y:180},{x:530,y:190},{x:530,y:370},{x:30,y:370}],
      sides: 4, color: 'rgba(230,126,34,0.2)', label: '4'
    }
  ];

  // Separatrices that can be removed during merging
  const separatrices = [
    { from: {x:250,y:160}, to: {x:270,y:30}, removedAt: 1 },  // between region 0 and 1
    { from: {x:250,y:160}, to: {x:320,y:180}, removedAt: 2 }, // between regions around center singularity
    { from: {x:320,y:180}, to: {x:400,y:30}, removedAt: 3 },  // between region 2 and 1
  ];

  // Singularities
  const sings = [
    {x: 250, y: 160, label: "+1"},
    {x: 320, y: 180, label: "+1"}
  ];

  // Merge state: which regions are merged
  const step = mergeStep;

  // After merging, redefine regions
  let visRegions = [];
  let visSeps = [];
  let visSings = [...sings];

  if (step === 0) {
    visRegions = regions;
    visSeps = separatrices;
  } else if (step === 1) {
    // Merge region 0 and 1 → new 4-sided region
    visRegions = [
      { pts: [{x:30,y:30},{x:400,y:30},{x:250,y:160},{x:30,y:170}], sides: 4, color: 'rgba(52,152,219,0.25)', label: '4' },
      regions[2], regions[3], regions[4], regions[5]
    ];
    visSeps = separatrices.filter(s => s.removedAt > 1);
  } else if (step === 2) {
    visRegions = [
      { pts: [{x:30,y:30},{x:400,y:30},{x:250,y:160},{x:30,y:170}], sides: 4, color: 'rgba(52,152,219,0.25)', label: '4' },
      regions[2],
      { pts: [{x:30,y:170},{x:250,y:160},{x:320,y:180},{x:280,y:260},{x:30,y:370}], sides: 4, color: 'rgba(241,196,15,0.25)', label: '4' },
      regions[5]
    ];
    visSeps = separatrices.filter(s => s.removedAt > 2);
  } else {
    visRegions = [
      { pts: [{x:30,y:30},{x:530,y:30},{x:530,y:190},{x:320,y:180},{x:250,y:160},{x:30,y:170}], sides: 4, color: 'rgba(52,152,219,0.2)', label: '4' },
      { pts: [{x:30,y:170},{x:250,y:160},{x:320,y:180},{x:530,y:190},{x:530,y:370},{x:30,y:370}], sides: 4, color: 'rgba(241,196,15,0.2)', label: '4' }
    ];
    visSeps = [];
    visSings = [];
  }

  // Draw regions
  for (const r of visRegions) {
    const d = 'M' + r.pts.map(p => `${p.x},${p.y}`).join(' L') + ' Z';
    svg.append("path").attr("d", d).attr("fill", r.color).attr("stroke", "#2c3e50").attr("stroke-width", 1.5);
    // Label
    const cx = r.pts.reduce((s,p) => s + p.x, 0) / r.pts.length;
    const cy = r.pts.reduce((s,p) => s + p.y, 0) / r.pts.length;
    svg.append("text").attr("x", cx).attr("y", cy).attr("text-anchor", "middle").attr("dominant-baseline", "middle")
      .attr("font-size", 18).attr("font-weight", "bold").attr("fill", "rgba(0,0,0,0.4)")
      .text(r.label + '-sided');
  }

  // Draw separatrices
  for (const s of visSeps) {
    svg.append("line").attr("x1",s.from.x).attr("y1",s.from.y).attr("x2",s.to.x).attr("y2",s.to.y)
      .attr("stroke","#e74c3c").attr("stroke-width",3).attr("stroke-dasharray","8,4");
  }

  // Draw singularities
  for (const s of visSings) {
    svg.append("circle").attr("cx",s.x).attr("cy",s.y).attr("r",7).attr("fill","#e74c3c");
    svg.append("text").attr("x",s.x).attr("y",s.y-14).attr("text-anchor","middle")
      .attr("font-size",11).attr("fill","#e74c3c").attr("font-weight","bold").text(s.label);
  }

  // Status text
  const statusTexts = [
    "Initial partition: 2 three-sided + 4 four-sided regions",
    "Merge step 1: removed separatrix → 3-sided region absorbed into 4-sided neighbor",
    "Merge step 2: removed another separatrix → second 3-sided region absorbed",
    "Final: all regions are 4-sided — ready for transfinite interpolation"
  ];
  svg.append("text").attr("x", W/2).attr("y", H - 8).attr("text-anchor", "middle")
    .attr("font-size", 13).attr("fill", "#333").text(statusTexts[step]);

  return svg.node();
}

Transfinite Bilinear Interpolation

Once every region is 4-sided, each is meshed independently using transfinite bilinear interpolation (TFI). Given a 4-sided region with boundary curves \(C_0(s)\) (bottom), \(C_1(s)\) (top), \(C_2(t)\) (left), \(C_3(t)\) (right) and corners \(P_{00}, P_{10}, P_{01}, P_{11}\), the interior mapping is:

\[ f(s,t) = (1-t)\,C_0(s) + t\,C_1(s) + (1-s)\,C_2(t) + s\,C_3(t) - \Big[(1-s)(1-t)P_{00} + s(1-t)P_{10} + (1-s)t\,P_{01} + st\,P_{11}\Big] \]

This is a Boolean sum of the boundary interpolants minus the bilinear correction at the corners. It maps the unit square \([0,1]^2\) to the curved quadrilateral, producing a structured grid of quadrilateral elements.

Drag the corners below to see how TFI generates a structured mesh inside an arbitrary 4-sided region:

8 x 8
Drag the orange corners to deform the block.
TipWhat to explore
  • Drag corners to create extreme quadrilateral shapes — notice how the scaled Jacobian drops when the quad becomes very distorted
  • A negative Jacobian (shown in red) means the mapping is inverted — the quad folds over itself, producing an invalid mesh element
  • Increase the grid density to see finer meshes — the quality metric stays the same regardless of resolution, since it depends only on the boundary shape

4. Results

4.1 Quality Comparison

The paper compares against unstructured quad meshes from GMSH on identical domains (~9000 quads each), using two metrics:

  • Minimum angle quality: the smallest angle in each quad (ideal = 90°)
  • Scaled Jacobian: measures element distortion (ideal = 1.0)

The block-structured approach achieves better quality on both metrics, primarily because of the low singularity count.

4.2 Examples

Symmetry preservation: The PDE-based approach naturally respects domain symmetry. A symmetric domain produces a symmetric partitioning — four 3-valent singularities and three 4-valent ones.

Airfoil: For a multi-element aerofoil (leading-edge flat, main body, trailing-edge flap), the method produces few singularities and generates boundary layers of quadrilaterals aligned orthogonally to the surface — exactly what CFD simulations need.

Multi-domain: When two domains share an interface curve, they can be partitioned separately (non-conforming) or jointly (conforming). Joint partitioning allows separatrices to cross between domains.

4.3 Computational Cost

“The time needed to partition the domain for the different examples presented in this paper varied from about one second, for most of them, to twenty seconds for the most complicated ones.”

The cost scales with the number of triangles in \(T_h\). The meshing step (transfinite interpolation) is negligible — even for millions of elements.

5. Connection to Our Work

This paper presents the same fundamental pipeline used in many frame field-based quad meshing methods:

Stage This Paper Ray et al. 2015
Field representation Representation vector \(\mathbf{u} = (\cos\theta, \sin\theta)\) Frame angle \(\theta\) with \(n\)-symmetry
Smoothness Dirichlet energy \(\int|\nabla\mathbf{u}|^2 dx\) Curvature \(\sum_e (\Delta\theta)^2\)
Solve method FEM + Lagrange multipliers for \(|\mathbf{u}|=1\) Direct linear solve in \((\cos 4\theta, \sin 4\theta)\)
Singularities Poincare index of \(\mathbf{u}\) Angle defect around triangles
Partitioning Separatrix tracing + transfinite interpolation MIQ parameterization + quad extraction

The key difference: this paper uses a representation vector \(\mathbf{u} \in \mathbb{R}^2\) with a unit norm constraint, while the Ray 2015 approach encodes cross field orientation through angle-based representations with 4-symmetry built in algebraically. Both arrive at the same cross field — just different numerical paths.

The NACA airfoil demo lets you explore several of these concepts interactively: boundary constraints, the linear solve, Jacobi smoothing, and streamline tracing through the resulting cross field.