Overview
In CAD-based turbomachinery design, blade profiles are defined by parametric curves (Bézier and B-Spline) whose control points the designer drags interactively. Geometric constraints — tangency to circles, lines, or other curves — must be maintained in real time as the user edits the shape, giving immediate visual feedback without manual re-solving.
This project developed a dual-solver algorithm that combines a novel linear matrix formulation with a nonlinear fallback optimizer. The linear solver uses minimum-norm control-point displacement to guarantee smooth transitions, while the nonlinear solver (SLSQP) handles over-constrained or degenerate cases that produce singular matrices.
The formulation was published at the Mechanical Engineering Research Symposium (MERS) 2025 and received the Best Poster Award. The production C++ implementation was integrated into a commercial turbomachinery design tool.
The Problem
Blade section editors in turbomachinery software let designers shape airfoil cross-sections using control points. These curves must satisfy geometric constraints — for example, a leading-edge curve must be tangent to both the pressure-side and suction-side curves as well as the stream-surface boundary. Moving any control point should automatically reposition the remaining points so that all constraints stay satisfied.
Two concrete use cases drove the work: (1) a diffuser airfoil editor where B-Spline curves must remain tangent to three circles (LE, TE, and thickness), and (2) a generic blade section editor where Bézier curves must be tangent to boundary lines and adjacent curves at their endpoints.
Initial Approach: Nonlinear Solver
Circle Tangency via Nonlinear Optimization
The first approach modeled the tangency constraints as a system of nonlinear equations. For a curve tangent to a circle, two conditions must hold simultaneously: the curve must pass through a point on the circle, and its tangent vector at that point must be perpendicular to the radius. With both the contact parameter and control-point positions as unknowns, this forms a nonlinear system.
A Python prototype using SciPy's gradient-free solver validated the concept for Bézier curves tangent to three circles. The approach was then generalized from 4 to arbitrary numbers of control points, and extended to B-Splines.
Generalized Algorithm (6+ Points)
The initial solver only handled 4 control points. A generalized version was built that works with any number of control points — essential for higher-degree curves that need more degrees of freedom.
B-Spline Extension
After Bézier validation, the algorithm was extended to B-Splines which use different basis functions and require a knot vector. This was needed for the diffuser airfoil editor where B-Splines define the profile.
Why a Nonlinear Solver Alone Isn't Enough
While the nonlinear optimizer found valid solutions, it had critical drawbacks for interactive use: control points would jump unpredictably between iterations, the solver would sometimes diverge or take too long for real-time feedback, and it offered no guarantee of minimal displacement — meaning the curve could change shape drastically even when a subtle adjustment would suffice.
For interactive operation where the solver runs on every mouse drag event, these issues made the pure nonlinear approach unsuitable. A fundamentally different formulation was needed.
Linear Solver: Minimum-Norm Formulation
Blade Section Line Tangency
The generic blade section editor required curves to be tangent to boundary lines and adjacent curves. This motivated a reformulation: instead of solving for absolute control-point positions, the algorithm solves for control-point displacements — the minimal movement needed to satisfy all constraints from the current configuration.
Expressing the problem in terms of displacements transforms it into a linear system, solvable via a minimum-norm solution that minimizes the total control-point movement. This guarantees smooth, predictable transitions during interactive editing.
Key Innovations
Directional tangency via cross-product constraint: Traditional tangent matching requires specifying a full tangent vector (direction + magnitude), consuming 2 DOF. By enforcing only the direction — via a cross-product that must equal zero — only 1 DOF is consumed. This was critical for avoiding over-constrained systems with few control points. The derivation uses a 2D rotation matrix to convert the cross-product into a linear dot-product form.
Start/end tangency as parametric lines: Rather than treating the second control point as an independent variable, it is re-parameterized as a point on a line through the first control point in the tangent direction. This introduces a scalar parameter instead of a 2D point, reducing the unknowns while maintaining the constraint structure as a linear system.
Iterative displacement for non-intersecting cases: When the curve doesn't yet intersect the target geometry, the algorithm samples candidate displacements and selects the one with minimum magnitude — then iterates until the constraint is met within tolerance. This bootstrapping step feeds the linear solver with the displacement vector it requires.
Displacement: No Intersection
When the curve doesn't intersect the target line, the algorithm samples points along the line, computes the displacement from each sample to the nearest curve point, and picks the minimum-displacement pair as the constraint input for the linear solver.
Displacement: With Intersection
When the curve already crosses the target, the algorithm identifies the longest overshoot on the wrong side and displaces that point back. This ensures convergence from both initial conditions — whether the curve starts inside or outside the constraint region.
Dual-Solver Architecture
Combining Linear and Nonlinear Solvers
The linear minimum-norm solver handles most interactive cases efficiently, but two situations require falling back to the nonlinear solver:
- Over-constrained systems: With only 4 control points and multiple tangency constraints, the linear coefficient matrix becomes singular. The nonlinear SLSQP optimizer can still find approximate solutions by treating the constraints as optimization targets.
- Circle tangency: Point-on-circle and tangent-to-circle constraints are inherently nonlinear (they involve squared radius terms). These are handled via equality constraints in the SLSQP formulation.
The architecture tries the linear solver first. If the coefficient matrix is singular or the error after solving exceeds tolerance, it automatically falls back to the nonlinear path. Both solvers share the same constraint specification API, making the fallback transparent.
Nonlinear Solver: SLSQP
The nonlinear path reformulates the problem as a constrained optimization: minimize total control-point displacement (or optionally control-polygon length for simpler shapes) subject to all geometric constraints as equality constraints.
Bounding boxes prevent the optimizer from exploring unreasonable regions of the solution space, and the initial guess is always the current control-point configuration — ensuring solutions stay close to the user's intent. The SLSQP algorithm (Sequential Least-Squares Quadratic Programming) was selected after testing several alternatives for its balance of speed and robustness on these small, dense problems.
Implementation
Software Architecture
The solver is implemented as a single C++ class (~3,100 lines) supporting Bézier, B-Spline, and NURBS curve types. The class exposes a clean constraint-specification API: callers add start/end tangency vectors, tangent-to-line, or tangent-to-circle constraints, then call InitSolver() for the initial solve and reCalculateControlPoints() on each interactive update.
Internally, the linear path constructs the coefficient matrix using Eigen for dense linear algebra — including Kronecker products, minimum-norm solutions via pseudoinverse, and column extraction for re-parameterized variables. The nonlinear path uses NLopt's SLSQP implementation with C-style callback wrappers and templated constraint/objective functions for type flexibility.
Key computational patterns include DOF analysis before solver selection, adaptive gradient vector updates when the user drags a specific control point, and post-solve validation that measures constraint error to decide whether to accept the result or retry with the fallback solver.
Python Blade Section Editor
All formulations were prototyped in Python first — SciPy for the nonlinear solver and NumPy for the matrix operations. An interactive blade section editor was built to validate real-time constraint maintenance during drag operations.
Two-Curve Constraint Demo
The final Python demo shows two curve segments with tangency constraints to circles, lines, and each other — validating that the combined solver handles the full set of constraints needed for production blade editing.
Publication
MERS 2025 — Best Poster Award
The linear formulation — specifically the directional tangency constraint derivation and the start/end tangency re-parameterization — was written up as a conference paper titled "Linear Matrix Formulation for Directional Tangency Constraints in Parametric Curve Design" and presented at the Mechanical Engineering Research Symposium (MERS) 2025 at the University of Moratuwa. The paper received the Best Poster Award.
The paper presents the mathematical framework in full detail; a separate page with the complete derivation is planned. The key contribution is showing that directional tangency — a constraint that traditionally requires nonlinear solving — can be expressed as a single linear row in the coefficient matrix using a 2D rotation matrix and cross-product identity.