Hao Hu

Mathias Kleiner

Jean-Philippe Pernot
email: jean-philippe.pernot@ensam.eu

Dr Mario Botsch

Dr Yongjie Jessica Zhang

Over-constraints detection and resolution in geometric equation systems ✩

Keywords:

This paper proposes an original decision-support approach to address over-constrained geometric configurations in Computer-Aided Design. It focuses particularly on the detection and resolution of redundant and conflicting constraints when deforming free-form surfaces made of NURBS patches. Based on a series of structural decompositions coupled with numerical analyses, the proposed approach handles both linear and non-linear constraints. The structural decompositions are particularly efficient because of the local support property of NURBS. Since the result of this detection process is not unique, several criteria are introduced to drive the designer in identifying which constraints should be removed to minimize the impact on his/her original design intent. Thus, even if the kernel of the algorithm works on equations and variables, the decision is taken by considering the user-specified geometric constraints. The method is illustrated on academic and industrial examples realized with our prototype software.

Introduction

Nowadays, designers rely on 3D CAD software to model sophisticated shapes based on free-form curves and surfaces. In industrial design, this geometric modeling step is often encapsulated in a larger Product Development Process (PDP) which may incorporate preliminary design, reverse engineering, simulation as well as manufacturing steps wherein several actors interact [START_REF] Falcidieno | Processing free form objects within a product development process framework[END_REF]. Actually, the final shape of a product often results from a long and tedious optimization process which tries to satisfy the requirements associated to the different steps and actors of the PDP. Requirements can be seen as constraints. They are generally expressed either with equations, a function to be minimized, and/or using procedures [START_REF] Gouaty | Variational geometric modeling with black box constraints and DAGs[END_REF]. The latter refers to the notion of black box constraints, not addressed in this paper, which focuses only on geometric constraints that can be expressed by linear or non-linear equations.

To satisfy the requirements, designers can act on variables associated to the different steps of the PDP. More specifically, in this paper, variables are supposed to be the parameters of the NURBS surfaces involved in the shape optimization process. To shape a free-form object defined by such surfaces, designers then have to specify the geometric constraints the object has to satisfy. For example, a patch has to go through a set of 3D points and satisfy to position constraints, the distance between two points located on a patch is fixed, two patches have to meet tangency constraints or higher-order continuity conditions, etc. Those geometric constraints give rise to a set of linear and non-linear equations linking the variables whose values have to be found. Due to the local support property of NURBS [START_REF] Piegl | The NURBS book[END_REF], the equations do not involve all the variables and some decompositions can be foreseen. Additionally, designers may express involuntarily several times the same requirements using different constraints thus leading to redundant equations. But, the designers may also involuntarily generate conflicting equations and may have to face over-constrained and unsatisfiable configurations.

Sometimes, over-constrained configurations can be solved by inserting extra degrees of freedom (DoFs) with the Boehm's knot insertion algorithm. As a consequence, many control points are added in areas where not so many DoFs are necessary [START_REF] Pernot | Fully free-form deformation features for aesthetic shape design[END_REF]. This uncontrolled increase of the DoFs impacts the overall quality of the final surfaces which become more difficult to manipulate than the initial ones. Furthermore, some structural over-constraints cannot disappear following this strategy and dedicated decision-support approaches have to be developed to identify and manage overconstrained configurations.

Unlike advanced 2D sketchers available in most commercial CAD software, and which can interactively identify the overconstraints during the sketching process, it is not yet completely possible to pre-analyze the status of 3D NURBS-based equation systems before submitting them to a solver. Thus, there is a need for developing a new approach for the detection and resolution of redundant and conflicting constraints in NURBSbased equation systems. This corresponds to the identification and treatment of over-constrained, well-constrained and underconstrained parts. In this paper, the treatment corresponds to the removal of constraints before solving. Once the constraints removed, the equation system often becomes under-constrained and the designer also has to add a requirement by mean of a function to be minimized so as to solve and find the values of the unknowns. This aspect is not part of the proposed approach but it will be discussed when introducing the results in which a particular functional is minimized.

Removing user-specified constraints is a perious step as the result do not fully satisfies what the designers have specified. Thus, it is not only important to develop an approach that is able to remove over-constraints, but also desirable to develop decisionsupport mechanisms which can help the designers identifying and removing the right constraints, i.e. the ones which preserve as much as possible the initial design intent.

This work contribution is to address these two difficult issues by proposing an original decision-support approach to manage over-constrained geometric configurations when deforming freeform surfaces. The algorithm handles linear as well as non-linear equations and exploits the local support property of NURBS. Based on a series of structural decompositions coupled with numerical analyses, the method detects and treats redundant as well as conflicting constraints. Since the result of this detection process is not unique, several criteria are introduced to drive the designer in identifying which constraints should be removed to minimize the impact on his/her original design intent. Thus, even if the kernel of the algorithm works on equations and variables, the decision is taken by considering the geometric constraints specified by the user at a high level.

The paper is organized as follows. Section 2 introduces the background and reviews the related works. Section 3 introduces the framework of our algorithm, details the principles and characteristics of its different steps and proposes criteria for evaluating its results. The proposed approach is then validated on both academic and industrial examples which are described in Section 4. Finally, Section 5 concludes this paper by discussing the main contributions as well as future work.

Background and related work

This section introduces how designers can specify their requirements within an optimization problem. It also analyses the existing methods used to detect structural or numerical over-constraints.

Modeling multiple requirements in an optimization problem

During the last decades, many deformation techniques have been proposed and it is not the purpose of this paper to detail all of them. Most of the time, when speaking of deformation techniques working on NURBS curves and surfaces, the goal is to find the position X of some control points so as to satisfy user-specified constraints which can be translated in a set of linear and/or nonlinear equations F (X ) = 0. Since the problem is often globally under-constrained, i.e. there are less equations than unknown variables, an objective function G(X ) also has to be minimized. As a consequence, the deformation of free-form shapes often results from the resolution of an optimization problem:

{ F (X ) = 0 min G(X ) (1)
For some particular applications, the optimization problem can also consider that the degrees, the knot sequences or the weights of the NURBS are unknown. However, in this paper, only the position of the control points are considered unknown. Depending on the approach, different objective functions can be adopted but they often look like an energy function which may rely on mechanical or physical models. The constraints toolbox can also contain more or less sophisticated constraints with more or less intuitive mechanisms to specify them.

Thinking to the PDP as well as to the needs for generating shapes which satisfy multiple requirements, one can notice that designers have access to three main parameters to specify their requirements and associated design intent within an optimization problem. They can effectively act on the unknowns X to decide which control points are fixed and which ones can move. In this way, they specify the parts of the initial shape which should not be affected by the deformation. Of course, designers can make use of the constraints toolbox to specify the equations F (X ) = 0 to be satisfied. Finally, designers can also specify some of their requirements through the function G(X ) to be minimized. For example, they can decide to preserve or not the original shape while minimizing an energy function characterizing the shape deformation.

However, most of the existing free-form shape deformation techniques do consider that the problem resulting from the set of equations F (X ) = 0 is under-constrained [START_REF] Elber | Geometric constraint solver using multivariate rational spline functions[END_REF][START_REF] Bartoň | Topologically guaranteed univariate solutions of underconstrained polynomial systems via no-loop and single-component tests[END_REF] and few attention has been paid to the analysis and processing of possible overconstraints. This paper proposes an approach to detect conflicting and redundant equations, and to help the designer in solving those issues by simply removing some constraints. However, Sections 3.4 and 4 discuss the possibility to fix more or less control points and thus modify the unknown vector X , as well as the possibility to modify the overall deformation behavior through the customization of the objective function G(X ) to be minimized.

Geometric over-constraints

Geometric over-constraints are classified into structural and numerical over-constraints [START_REF] Sridhar | Algorithms for the structural diagnosis and decomposition of sparse, underconstrained design systems[END_REF]. Structural over-constraints can be detected from an analysis of the DoFs, at the level of either the geometry or the equations. Numerical over-constraints are usually determined from an analysis of the solvability of the equations system. Since our approach is based on equations, both aspects are to be defined.

Structural over-constraints

Jermann et al. give a general definition of structurally overconstrained, well-constrained and under-constrained equation systems at a rather macro level and considering the dimension of the space [START_REF] Jermann | Decomposition of geometric constraint systems: a survey[END_REF]. This definition has been here adapted to system of equations where the system is expected to be fixed with respect to a global coordinate system. Definition 1. The degree of freedom DoF (v) of a geometric entity v is the number of independent parameters that must be set to determine its position and orientation. For example, in 2D space, it is equal to 2 for points and lines. For a geometric constraints system G with a set V of geometries, the degree of freedom of all the geometries is DoFs = ∑

v∈V DoF (v). Finally, the above DoF-based counting definitions compare the number of equations to the number of variables of a system. However, they do not cover cases such as geometric redundancies induced by geometric theorems. In order to cover those situations, algebraic definitions are introduced.

Numerical over-constraints

The previous structural definitions cannot distinguish redundant and conflicting constraints. However, from the algebraic point of view, it is ascertain that this could be handled properly by Grobner basis or Wu-Ritt methods [START_REF] Chou | Ritt-Wu's decomposition algorithm and geometry theorem proving[END_REF]. These methods are commonly used in abstract algebra and requires strong mathematical fundamentals to understand. Definition 5. Let G = (E, V , P) be a geometric constraints system, where E is a set of equations, V is a set of variables and P is a set of parameters. The set of solutions to G is denoted Sol(G). A geometric constraints system is inconsistent iff Sol(G) = ∅ and is consistent iff Sol(G)̸ = ∅. Definition 6. Let G = (E, V , P) be a consistent geometric constraints system. Let G ′ = (E ∪ E c , V , P ′ ) be an inconsistent geometric constraints system, where E c is a set of equations forming a constraint C = {E c | E c ∩ E = ∅} and P ⊂ P ′ . As a result, C is a conflicting constraint with respect to G. Definition 7. Let G = (E, V , P) be a consistent geometric constraints system. Let G ′ = (E ∪ E r , V , P ′ ) be a consistent geometric constraints system, where E r is a set of equations forming a constraint R = {E r | E r ∩ E = ∅} and P ⊂ P ′ , and Sol(G) is the same as Sol(G'). As a result, R is a redundant constraint with respect to G. Definition 8. Let G = (E, V , P) be a weakly connected geometric constraints system (its components are weakly connected [START_REF] Diestel | Graph theory[END_REF]) which can be decomposed into the following two subsystems: It should be noted that the set of basis constraints and numerical over-constraints of a given system is not unique. Thus, decisionsupport mechanisms and criteria have to be defined to help designers identifying the right redundant and conflicting constraints to be removed.

G b = (E b , V , P) and G o = (E o , V , P) with {E = E b ∪E o , E b ∩E o = ∅}. If any constraint E oi in E o is

Modeling geometric constraints system

As previously discussed, a geometric constraints system can be described either at the level of the equations or at the level of the geometry. On the one hand, for a system of equations, there exists algebraic methods able either to directly address consistency problems [START_REF] Cox | Ideals, varieties, and algorithms: An introduction to computational algebraic geometry and commutative algebra[END_REF] or to analyze the structure indirectly based on a bipartite graph where two classes of nodes represent variables and equations independently [START_REF] Bunus | A debugging scheme for declarative equation based modeling languages[END_REF]. On the other hand, for modeling at the level of the geometry, two types of graph are mostly used: either bipartite graphs with two classes of nodes representing geometric entities and constraints separately [START_REF] Hoffmann | Geometric constraint decomposition[END_REF], or constraint graphs with nodes representing the geometric entities and edges representing constraints [START_REF] Gao | Solving geometric constraint systems. ii. a symbolic approach and decision of rc-constructibility[END_REF][START_REF] Hoffman | Decomposition plans for geometric constraint systems, part i: performance measures for cad[END_REF]. However, considering NURBS-based constraint systems is not straightforward since there exists several types of variables contributing to the shape deformation. On one hand, system of equations enable a sound way of modeling where parameters like knots and weights can be set as variables. On the other hand, graph modeling at the level of the geometry is currently limited to variables like the coordinates of the control points and their associated weights. Representing variables such as degrees, knots, values of u and v parameters using constraint graphs at the level of the geometry is not convincingly demonstrated in the literature. In the work of Lesage [START_REF] Lesage | Un modèle dynamique de spécifications d'ingénierie basé sur une approche de géométrie variationnelle[END_REF], vectors between control points are used to represent NURBS objects as well as the geometric constraints such as incidence or tangency. Nevertheless, his method is restricted to cases where control points are only unknowns and is not general enough comparing to equation-based modeling.

According to the previous definitions, detection methods can be classified in two categories [START_REF] Sridhar | Algorithms for the structural diagnosis and decomposition of sparse, underconstrained design systems[END_REF]: structural over-constraints and numerical over-constraints detection.

Structural over-constraints detection

Many structural over-constraints can be identified while counting DoFs during the system decomposition process. Among the decomposition methods, those that are helpful for finding structurally over-constrained subparts have been identified and classified in three categories according to the type of subsystems.

Specific patterns

This category is based on the fact that, when considering their engineering drawings, most systems can be decomposed while recognizing specific patterns constructed by ruler and compass. Recursive division is first proposed by Owen to handle 2D constraint systems where only distance and angle constraints are involved [START_REF] Owen | Algebraic solution for geometry from dimensional constraints[END_REF]. He introduced several rules for redundancy detection, including rules to detect over-rigid subparts and rules to check angle redundancies [START_REF] Owen | Constraints on simple geometry in two and three dimensions[END_REF]. Fudos and Hoffman adopted the bottom-up graph reduction method that works well with over-and well-constrained systems in two dimensions [START_REF] Fudos | A graph-constructive approach to solving systems of geometric constraints[END_REF]. These methods are polynomial in time but not general enough due to the limited repertoire of patterns, which cannot cover all types of geometric configurations.

Structural rigidity

This class gathers together methods that decompose a system into structural rigid subsystems. Methods vary depending on the adopted structural-rigidity definition as well as on the corresponding search algorithms. By modifying incremental network maximum flow theory, Hoffmann et al. developed the Dense algorithm to identify 1-well-constrained subgraph [START_REF] Hoffmann | Making constraint solvers more usable: overconstraint problem[END_REF]. Their definition of structural-rigid subsystem derives from Laman's theorem on the characterization of the rigidity of bar frameworks [START_REF] Graver | Combinatorial rigidity[END_REF]. However, the definition does not handle properly constraints such as incidences and parallelisms, which are widely used in CAD systems. Jermann et al. modified their definition by introducing the notion of degree of rigidity (DoR) to replace the dimension D [START_REF] Jermann | Algorithms for identifying rigid subsystems in geometric constraint systems[END_REF]. The difference between the two lies on the fact that the value of DoR varies with subsystems while D remains constant whatever the subsystems are (e.g. in 2D, D = 3 and in 3D, D = 6). Based on the same network flow theory, his algorithm Over-rigid deals properly with specified parallel/incidence constraints but is still restricted to generic conditions where incidence degeneracies due to geometric theorems like co-linearities, co-planarities are forbidden [START_REF] Jermann | Algorithms for identifying rigid subsystems in geometric constraint systems[END_REF].

Maximum matching

Those methods recognize structural over-constrained patterns by directly comparing the DoFs and DoCs of a system (or subsystem) without considering the dimension-dependent constant D. The Dulmage-Mendelsohn (D-M) decomposition algorithm allows for decomposing a system of equations into over-constrained, well-constrained and under-constrained subsystems [START_REF] Dulmage | Coverings of bipartite graphs[END_REF]. It has been used for debugging in equation-based modeling systems such as Modelica [START_REF] Bunus | A debugging scheme for declarative equation based modeling languages[END_REF]. Additionally, it computes a directed acyclic graph (DAG) which provides a solving order among the strongly connected components (SCC) of the system. Serrano has been interested in using graph-theoretic algorithm to prevent overconstrained systems where all constraints and geometric entities are of DoF one [START_REF] Serrano | Constraint management in conceptual design[END_REF]. Latham et al. extended the work of Serrano by proposing maximum b-weighted matching to identify overconstraints with arbitrary DoFs [START_REF] Latham | Connectivity analysis: a tool for processing geometric constraints[END_REF].

A structurally over-constrained part contains either redundant or conflicting constraints. But the remaining parts may also contain numerical over-constraints [START_REF] Podgorelec | Dealing with redundancy and inconsistency in constructive geometric constraint solving[END_REF]. This is because structure-based analysis ignores numerical information of a system. Hence, to further identify subtle over-constraints like geometric redundancy, numerical methods have to be adopted.

Numerical over-constraints detection

Any geometric constraint can be transformed into a set of algebraic equations [START_REF] Hoffmann | Geometric constraint decomposition[END_REF]. Therefore, geometric over-constraints are equivalent to a set of inconsistent (i.e. conflicting) or redundant equations. Here, numerical detection methods have been classified into two categories according to the type of constraints.

Linear over-constraints detection

Gauss elimination, LU factorization with partial pivoting and QR factorization with column pivoting have been successfully adopted to find redundant/conflicting equations as well as spanning group in systems of linear equations [START_REF] Strang | Linear algebra and its applications[END_REF]. Light and Gossard applied Gauss elimination to compute the rank as well as to further identify invalid equations [START_REF] Light | Variational geometry: a new method for modifying part geometry for finite element analysis[END_REF]. Serrano extended their work to check existence of over-constraints within strongly connected components of a system of equations [START_REF] Serrano | Constraint management in conceptual design[END_REF]. These methods enable stable and fast detections but are limited to linear cases.

Non-linear over-constraints detection

Symbolic methods are theoretically reliable but the time complexity is exponential. Kondo used Grobner basis method to test dependency among 2D dimension constraints [START_REF] Kondo | Algebraic method for manipulation of dimensional relationships in geometric models[END_REF]. Gao and Chou introduced Wu-Ritt's decomposition algorithm to determine whether a system is over-constrained [START_REF] Gao | Solving geometric constraint systems. ii. a symbolic approach and decision of rc-constructibility[END_REF]. However, both methods do not directly find the spanning groups of over-constraints.

Optimization methods have been used to address constraints satisfaction problems, which works well for under-constrained systems [START_REF] Ge | Geometric constraint satisfaction using optimization methods[END_REF].

Jacobian matrix analysis methods enable a faster detection by studying Jacobian structure of equations. However, they are not able to distinguish redundant and conflicting constraints. The main difference between these methods is the configuration where the Jacobian matrix should expand. If the system is solvable, Haug proposed to perturb the common root and recalculate the rank once the Jacobian matrix is rank deficient [START_REF] Haug | Computer aided kinematics and dynamics of mechanical systems[END_REF]. However, if the system is non-solvable, Foufou et al. suggest a Numerical Probabilistic Method (NPM), which analyzes the Jacobian matrix at random configurations [START_REF] Foufou | Interrogating witnesses for geometric constraint solving[END_REF]. However, there is a risk that the Jacobian matrix is rank deficient at the chosen points but is full rank everywhere else. Therefore, NPM is convenient in computation but may lead to incorrectly detected over-constraints. Instead of randomly selecting configurations, Michelucci et al. suggested to study the Jacobian structure at witness configurations where incidence constraints are satisfied [START_REF] Michelucci | Geometric constraint solving: the witness configuration method[END_REF]. The witness configuration and the target configuration share the same Jacobian structure. As a consequence, all the over-constraints are identified. More recently, Moinet et al. developed tools to identify conflicting constraints through analyzing the witness of a linearized system of equations [START_REF] Moinet | Defining tools to address over-constrained geometric problems in computer aided design[END_REF]. Their approach has been applied to the well-known Double-Banana test case on which our approach will also be tested in Section 4.

From the above discussion, it is clear that different methods are capable of handling certain geometric constraints systems. However, no method is able to cover all cases perfectly within our criteria. Thus, in this paper, a new approach which couple structural as well as numerical methods is proposed.

A generic approach coupling structural decompositions and numerical analyses

This section describes our approach for detecting and treating redundant and conflicting geometric constraints. The main idea is to decompose the system of equations into smaller blocks that can be analyzed iteratively using dedicated numerical methods. The overall framework and algorithm are first introduced before detailing the different steps involved.

Overall detection framework

The overall framework has been modeled in Fig. 1. It is based on three nested loops: the structural decomposition into connected components (CC ); the structural decomposition of a CC into its subparts (G1, G2, G3) and its corresponding DAG of strongly connected components (SCC ); the iterative numerical analysis of these SCC . Pseudo-code for the main procedures is provided in Section 3.2.

Loop among connected components

The system of equations (SE) is initially represented by a graph structure G, where nodes correspond to variables and edges to equations. The structure is first decomposed into n connected components {CC 1 , . . . , CC n } using Breadth First Search (BFS) [START_REF] Leiserson | A work-efficient parallel breadth-first search algorithm (or how to cope with the nondeterminism of reducers)[END_REF]. Such a decomposition is made possible thanks to the local support property of NURBS or simply when using constraints which decouple what happens along the x, y and z directions of the reference frame (e.g. position or coincidence constraints). As a result, geometric over-constraints can be detected separately for each CC i .

Loop among subparts obtained by D-M decomposition

The D-M decomposition is used to structurally decompose CC i into a maximum of three subparts: G i1 (over-constrained subpart), G i2 (well-constrained subpart) and G i3 (under-constrained subpart). Each subpart (if it exists) will be analyzed iteratively using the third nested loop explained below.

However, a single pass of the third loop on each G ij is not sufficient. Indeed, any pass may lead to the removal of constraints, which modifies the structure of the CC i and thus requires to apply D-M decomposition again after the pass to obtain updated subparts. The superscript d is used to note that CC d i (resp. G d ij ) refers to CC i (resp. G ij ) after its dth D-M decomposition. Although the number of passes required is unknown in advance, it is guaranteed that the process will converge to a state where only one subpart G i3 is left. In other words, constraints will be either removed or moved to the third subpart along the process. For each red block, and as shown in the top-left part of the Fig. 1, an appropriate numerical method (numFindRC in the figure and the pseudo-code) tries to find redundant (R) or conflicting (C) constraints. These over-constraints are then removed from the currently analyzed connected component CC d ij . If the block is solvable, its (unique) solution is propagated to dependent blocks, which may lead to additional redundant or conflicting constraints being detected and removed from CC d ij . Once all red blocks have been analyzed, this part of the DAG (potentially turning blue blocks into red ones) is recomputed until all blocks are analyzed. However, it is not needed to recompute the D-M decomposition on the whole CC d ij , it is sufficient to recompute only the maximum matching for the current subpart by calling linkedSCC again. The superscript m is used to note that G dm ij (resp. SCC dm ijk ) refers to G d ij (resp. SCC d ijk ) after its mth matching. Although the number of passes required is unknown in advance, it is guaranteed that the process will converge to a state where there are either no blocks left, or these blocks only contain variables (and that is only possible for the third subpart G d i3 ). In other words, constraints and variables are removed until we obtain an under-constrained system with multiple solutions, meaning no more propagation is possible. In that last step, as shown in the bottom-left part of the figure, the remaining system is analyzed for numerical conflicts and proceeds with the next connected component.

Pseudo-code

This section provides the pseudo-code for the two main procedures of the approach, surrounded by dotted rectangles on Fig. 1.

Algorithm 1 Structural decomposition

1: SE ← System of Equations 2: G ← Graph(SE) 3: [CC 1 , • • • , CC n ] ←BFS(G) 4: for i = 1 to n do 5: [G 1 i1 , G 1 i2 , G 1 i3 ] ←DM(CC i ) 6: CC 1 i ← CC i 7:
for j = 1 to 3 do

8:

d ← 1 9: continue←True 10:

while continue and G d ij ̸ = ∅ do 11: continue,CC d+1 i ← findRC(CC d i , G d ij )
12:

d ← d + 1 13: [G d i1 , G d i2 , G d i3 ] ←DM(CC d i )
14:

end while

15:

end for 16: end for

17: return [CC d 1 , • • • , CC d n ]

Strongly connected components analysis

This section discusses the techniques used to analyze the strongly connected components SCC dm ijk . This corresponds to numFindRC function of Algorithm 2 (Section 3.2) used to find redundant (R) and conflicting (C) constraints of a component if they exist. Otherwise, the component is solved and the solutions are propagated to the whole system.

Depending on the type of constraints, i.e. either linear or nonlinear ones, methods differ and are presented in the next subsections. The following notation A[i : j, l : k] is used to define the matrix obtained by slicing the ith to jth rows, and the lth to kth columns of A.

CC d i ← removeRCfromCC(CC d i , [R, C ])

Linear system

In the proposed approach, the QR factorization with column pivoting is used to detect linear over-constraints. QR factorization with an optional column permutation P, triggered by the presence of a third output argument, is useful for detecting singularity or rank deficiency. Fig. 2 shows the overall detection process. The horizontal colored straight lines correspond to the linear equations of the system A.x = b to be solved, where A has a dimension m × n. Here, it is assumed that the rank of the system is r, which means that there are r independent equations with m > r. The rank is computed using SVD, which is relatively stable compared to other methods.

As for QR factorization, columns are exchanged at the start of the kth stage to ensure that

  A (k) k (k : m)    2 = max j⩾k   A (k) j (k : m)    2 (2)
where

A (k) j (k : m) = A[k : m, j].
At each step of the factorization, the column of the remaining un-factored matrix with largest norm is used as the basis for that step and is moved to the lead position [START_REF] Golub | Matrix computations[END_REF]. This ensures that the diagonal elements of R occur in decreasing order and that any linear dependence among the columns is certainly revealed by examining these elements. Permutation matrix P rearranges the columns of A t so that the columns appear in the decreasing order of their norm.

The first r columns of A t .P are the basis constraints of A t and the first r columns of Q form an orthogonal basis (Fig. 2b). Since the remaining mr columns are linearly dependent on the first r columns [START_REF] Dongarra | Proceedings of the fourth SIAM conference on parallel processing for scientific computing[END_REF], they are the over-constraints. The rank r also corresponds to the number of non-zero values of diagonal elements of R.

To find linear dependencies between the columns, the following deduction is needed. First, the matrix Q (:, 1 : r) is inverted using the following equation: and is then used in the following equation:

A t (:, 1 : r) = Q (:, 1 : r).R(1 : r, 1 : r)
A t (:, r + 1 : n) = Q (:, 1 : r).R(1 : r, r + 1 : n)
thus providing the following relationship between the two sliced matrices A t (:, r + 1 : n) and A t (:, 1 : r):

A t (:, r + 1 : n) = A t (:, 1 : r).R(1 : r, 1 : r) -1 R(1 : r, r + 1 : n)
Finally, to identify the redundant and conflicting equations, the new b vector after factorization is redefined as follows:

b new = b(r + 1 : n) -b(1 : r).R(1 : r, 1 : r) -1 R(1 : r, r + 1 : n)
Redundant and conflicting equations are further distinguished by comparing the value of the last mr elements of b new with 0.

Non-linear system

When considering a system of non-linear equations, a twophases identification process is used. First, the Witness Configuration Method [START_REF] Michelucci | Geometric constraint solving: the witness configuration method[END_REF] is used to find all the over-constraints (phase I), and Grobner Basis or Incremental Solving is then applied to further distinguish redundant and conflicting constraints (phase II).

Phase I. Taking advantage of the method proposed by Moinet et al. [START_REF] Moinet | Defining tools to address over-constrained geometric problems in computer aided design[END_REF], a generic witness configuration is generated from the initial shape of the object to be deformed (step 1). Effectively, in our case, the variables x are the positions of the control points which have an initial location x (0) before deformation. Then, QR factorization with column pivoting is used to analyze this witness configuration (step 2). As a result, the sequence of equations is reordered. The first r (rank of the Jacobian matrix) equations are independent while the remaining ones are the over-constraints. In Fig. 3, the curved colored lines represent non-linear equations.

Phase II. To further distinguish redundant and conflicting constraints, either Grobner Basis or Incremental Solving are applied. In our algorithm, this choice relies on the number of equations. If the number of equations m ≤ 10, Grobner Basis is preferred [START_REF] Lamure | Qualitative study of geometric constraints[END_REF]. Oth- erwise, Incremental Solving is chosen. To explain the two methods, let us assume that the following constraints system is available after Phase I (Fig. 4): where the Eq. ( 4) to r are the basis constraints and the equations (r + 1) to m are the over-constraints.

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ f 1 (x 1 , x 2 , . . . , x n ) = 0 f 2 (x 1 , x 2 , . . . , x n ) = 0 . . . f r (x 1 , x 2 , . . . , x n ) = 0 f r+1 (x 1 , x 2 , . . . , x n ) = 0 . . . f m (x 1 , x 2 , . . . , x n ) = 0 (3)
The Incremental Solving method incrementally insert the overconstraint f j = 0, j ∈ {r + 1, . . . , m}, in the set of basis constraints thus forming a new group of equations {f 1 = 0, . . . , f r = 0, f j = 0}. If the new group is solvable, then the equation f j = 0 is redundant, otherwise it is conflicting. Of course, the basis constraints are always solvable. When Grobner basis [START_REF] Cox | Ideals, varieties, and algorithms: An introduction to computational algebraic geometry and commutative algebra[END_REF] are used, the method first computes the reduced Grobner basis rgb r of the ideal ⟨f 1 , . . . , f r ⟩. Since the set of equations are solvable, rgb r ̸ = {1}. Then, a loop on all the over-constraints f j = 0, j ∈ {r + 1, . . . , m}, starts and for each over-constraint the reduced Grobner basis rgb r+j of the ideal

⟨f 1 , . . . , f r , f j ⟩ is computed. If rgb r+j ≡ rgb r , then f j = 0 is a redundant equation. If rgb r ⊂ rgb r+j , then f j = 0 is an independent equation. Finally, if rgb r+j = {1}, then f j = 0 is a conflicting equation.

Validation and evaluation of the solutions

Section 2.1 has introduced the multiple ways to model requirements within an optimization problem by specifying an unknown vector X , the constraints to be satisfied F (X ) = 0 and the function G(X ) to minimize.

The approach described in this section allows for the identification of redundant and conflicting equations. Correctness is ensured since it consists of a fixed-point algorithm that only stops when the system is solvable. Additionally, any removed equation is guaranteed to be either conflicting or redundant with the remaining set. It has thus been shown that the set of equations F (X ) = 0 can be decomposed in two subsets: F b (X ) = 0 containing the basis equations, and F o (X ) = 0 the over-constrained ones.

To stay close to the requirements the designer has in mind, the proposed approach then moves from the equations level to the constraints level. Thus, the geometric constraints associated to the equations F o (X ) = 0 are analyzed and all the equations related to those constraints are gathered together in a new set of equations Fo (X ) = 0. Of course, the equations F o (X ) = 0 are included in the set of equations Fo (X ) = 0. Finally, the equations related to constraints which are neither conflicting nor redundant form the other set Fb (X ) = 0. This transformation allows working at the level of the constraints and not at the level of the equations. This is much more convenient for the end-user interested in working at the level of geometric requirements.

Since this decomposition is not unique, it gives birth to various potential final solutions (interactive decomposition is out of scope of this paper). Therefore, several criteria are now introduced to evaluate these solutions according to the initial design intent. To be able to characterize the quality of the obtained solutions, the set of user-specified parameters P is introduced. This set gathers together all the parameters the designer can introduce to define the constraints his/her shape has to satisfy. For example, the distance d imposed between two points of a NURBS surface is a parameter characterizing a part of the design intent. Then, the idea is to evaluate how much the solutions deviate from the initial design intent and notably in terms of the parameters P.

To do so, the optimization problem containing the basis constraints is solved:

{ Fb (X ) = 0 min G(X ) (4)
and the solution X ′ is then used to evaluate the unsatisfied overconstraints Fo (X ′ ) as well as the real values P ′ of the user-specified parameters P. For example, if the user-specified distance d between the two patches cannot be met, then the real distance d ′ will be measured on the obtained solution. From this solution, it is then possible to evaluate three quality criteria:

• Deviation in terms of parameters/constraints: this criterion aims at measuring how far/close the real values P ′ of the parameters are from the user-specified parameters P. This criterion helps understanding if the design intent is preserved in terms of parameters and consequently in terms of constraints.

df = ∑ i ⏐ ⏐ P ′ i -P i ⏐ ⏐ ∑ i |P i | (5)
• Deviation in terms of function to minimize: this criterion directly evaluates how much the objective function G has been minimized. Here, the function is simply computed from the solution X ′ of the optimization problem. To preserve the design intent this value is to be minimized. Thus, it can be used to compare the solutions between them.

dg = G(X ′ ) (6)
• Degree of near-dependency: rank deficiency of the Jacobian matrix at the witness clearly reveals the dependencies between constraints. However, for NURBS-based equation systems, the constraints can be independent but near to be dependent. In this case, the Jacobian matrix of Fb (X ) at the solution point X ′ is ill-conditioned and the corresponding solution can be of bad quality. The third criterion thus evaluates the condition number (cond) of the Jacobian matrix as a measure of near-dependency [START_REF] Kincaid | Numerical analysis: mathematics of scientific computing[END_REF]:

cond = cond(J JF b (X ′ )) (7)
Fig. 5. Initial geometry of the Double-Banana as described in [START_REF] Moinet | Defining tools to address over-constrained geometric problems in computer aided design[END_REF].

Finally, even if those criteria characterize the quality of the solution X ′ with respect to the design intent, they have not been combined in a unique indicator. Thus, the results of the next section will be evaluated by analyzing and comparing those three criteria for each solution.

Results and discussion

This section presents two configurations on which the proposed over-constraints detection and resolution technique has been tested. The first one concerns the academic Double-Banana testing case widely studied in the literature. It has been used to be able to compare our solution to the ones generated by others. The second example is more industrial and concerns the shaping of a glass composed of several NURBS patches.

Double-Banana testing case

The variables X , the constraints F (X ) = 0, and the parameters P of the Double-Banana testing case are exactly the same as the ones tested by Moinet et al. [START_REF] Moinet | Defining tools to address over-constrained geometric problems in computer aided design[END_REF]. The only difference is that they are using a coordinate-free formulation whereas ours is cartesianbased. Here, the objective is to find the position of the 8 nodes of a 3D structure so that the length of the 18 edges satisfy userspecified dimensions. Fig. 5 illustrates the Double-Banana in its initial configuration.

The Double-Banana configuration contains only one connected component as revealed by BFS. The structural analysis using D-M decomposition shows that it is under-constrained and our algorithm then follows the bottom part of the Fig. 1. WCM analysis is used within our numFindRC function and an over-constraint is detected. More specifically, the equation e9 is here detected. Using our Incremental Solving approach, the equation is further characterized as a conflicting one. The equation e9 is therefore removed and the system is solved using the initial position of the nodes as initial values of the variables. Using the results, the equation e9 is then reevaluated and the associated parameter is compared to the user-specified value. In the present case, e9 is not satisfied since it is equal to 44.47 compared to the initial user-specified requirement of 45. Thus, the deviation from the design intent is df = 0.53/45.

Our algorithm gives a solution that is much closer to the initial design intent than the algorithm of Moinet et al., and the remaining system is less ill-conditioned after removing the conflicting constraint (Table 1). Actually, the algorithm of Moinet et al. identifies e18 as a conflicting equation and its removal induces a deviation df = 4.38/32 from the initial design intent.

Sketching a 3D glass

In this example, the idea is to show how the proposed over-constraints detection and resolution approach can support the sketching of a 3D glass composed of 4 connected NURBS patches. The designer sketches his/her design intent and associated requirements. Here, the objective is to modify the upper part of the glass by specifying the following elements:

• Variables: Each patch has a degree 5 × 5 and has a control polygon made of 16 × 6 control points whose coordinates are the variables of our optimization process (Fig. 6a). Since the objective is to modify the upper part of the glass, the designer selects how many rows of control points are to be blocked and how many can move. For example, if the designer wishes to free the upper row of control points of the four patches, then there will be 6 × 4 × 3 = 72 variables in the unknown vector X . The results will be illustrated with 4 and 5 rows free to move.

• Constraints: Three types of constraints are used to specify how the shape of the 3D glass has to evolve:

-Position: 4 position constraints are added to the four end points of the patches along the upper boundary curves. As shown in Fig. 6c, the green points of the patches need to move to new positions in the 3D space. They are labeled from 1 to 4 and they generate 4 × 3 = 12 linear equations labeled from 1 to 12 (Table 2). -Distance: 2 distance constraints are defined between the opposite sides of the patches (Fig. 6d). They are labeled 5 and 6 and they generate 2 × 1 non-linear equations labeled 13 and 14 (Table 2). -Coincidence and tangency: 8 coincidence and 8 tangency constraints are specified to maintain the continuity between the upper parts of the patches during the deformation (Fig. 6b). They are labeled from 7 to 22 and they generate 8 × 3 linear and 8 × 3 non-linear equations labeled from 15 to 62 (Table 2). The nonlinearity comes from the use of the vector product to express the collinearity of normals.

Overall, there are 22 geometric constraints generating 62 equations in the set F (X ) = 0. Some of those constraints are conflicting and it is the purpose of this section to try to see how our algorithm can detect and remove them without affecting too much the design intent.

1 0 0 1 1 1 1 2 1 0 0 1 1 1 3 1 0 1 1 1 0 4 0 1 1 0 1 1 5 0 1 1 1 0 1 6 1 1 0 0 1 1 7 1 1 1 1 0 0 8 1 1 1 0 1 0 9 1 1 0 1 0 1
Fig. 6. Initial sketch of glass geometry.

• Objective function: Since the proposed approach removes the identified over-constraints, the resulting system of equations Fb (X ) = 0 (Section 3.4) may become under- constrained and a function G(X ) has to be minimized. Here, the idea is to make use of the approach of Pernot et al. to define two types of deformation behavior [START_REF] Pernot | Fully free-form deformation features for aesthetic shape design[END_REF]: either a minimization of the variation of the shape (G 1 (X )) between the initial and final configurations, or minimization of the area of the final shape (G 2 (X )). In terms of design intent, the first one tends to preserve the initial shape of the glass, whereas the second forgets the initial shape and tends to generate surfaces similar to tensile structures.

As revealed by BFS, the system can be decomposed into 2 + 24 × N rows connected components CC i where N rows is the number of rows free to move. Among them, only two components CC 1 and CC 2 contain both variables and equations while the others contain only variables (Table 2). The analysis of those two components gives rise to the identification of 2 conflicting equations which correspond to either the position or distance constraints. Since the result of the detection process is not unique, 9 configurations are obtained and are gathered together in Table 3. Here, one has to remember that even if the detection process identifies conflicting equations, our algorithm removes the constraints associated to those equations. For example, configuration 1 considers that the two distance constraints (one between patches P1 and P4, and the other between P2 and P3) are to be removed (0 in the table) and the 4 position constraints are kept (1 in the table).

All configurations are then solved while acting on both the number of upper rows to be fixed (N rows = 4 or 5), and the objective function to be minimized (either G 1 (X ) or G 2 (X )). The results are gathered together in Tables 4 and5. Each configuration is evaluated through the three previously introduced criteria dg, df and cond. Some solutions are shown in Fig. 7.

One can first notice that depending on the configurations, the deviation df on the constraints varies. For example, with N rows = 4 and while minimizing G 1 (X ), the configuration 7 generates a solution that is closer to the design intent than configuration 6 (0.10684 < 0.12607 in Table 4). For configuration 3, it is clear that the deviation to the design intent in terms of constraints is more important when minimizing the area of the final surface than when minimizing the shape variation (0.2288 > 0.10179 in Table 4). This is clearly visible on Fig. 7c 1 and c 2 .

But the deviation dg i on the objective function to be minimized also varies. While considering the minimization of the shape variation, one can see that configuration 3 is less interesting than configuration 1 in the sense that it minimizes less the shape variation (15459.52 > 13801.04 in Table 4).

Finally, for a given configuration, one can notice that when the number of free rows increases, i.e. when there is more freedom, the objective function decreases and the solution is therefore closer to the design intent. This is visible when comparing values from Tables 4 and5. Thus, the selection of the variables X are also important when setting up the optimization problem.

Conclusion and future work

In this work, an approach for finding all over-constraints in free-form geometric configurations has been introduced. It relies on a coupling between structural decompositions and numerical analysis. The process and its algorithm have been described and analyzed with results on both academic and industrial examples. The approach has several benefits: it is able to distinguish between redundant and conflicting constraints; it is applicable on both linear and non-linear constraints; and it applies numerical methods on small sub-blocks of the original system, thus allowing to scale to some large configurations. Additionally, since the set of overconstraints of a system is not unique, it has been shown that our approach is able to provide different sets depending on the selected structural decomposition, and proposed criteria to compare them and assist the user in choosing the constraints he/she wants to remove. Even if the kernel of the algorithm works on equations and variables, the decision is taken by considering the geometric constraints specified by the designer at a high level.

A number of perspectives stem from this work. First, an automatization of the process should assist the designer in selecting the set of over-constraints that less deviate from his/her original design intent. As it is, the designer has access to three main criteria (dg, df , cond) which can be difficult to analyze for a non-expert. Thus, higher-level criteria should be imagined on top of those ones. Second, the approach can be made interactive, i.e. allowing the designer to select between the different conflicting sets along the process, or even modify the faulty constraints. Finally, it is planned to extend this work so that it can be used to detect and explain geometric configurations which, even when solvable, result in poor quality designs.

Fig. 1 .

 1 Fig. 1. Overall framework composed of three nested loops defining the main structure of the detection algorithm.

3. 1 . 3 .

 13 Loop among strongly connected components In addition to the subparts, D-M decomposition also provides a DAG for each CC d i . Nodes of this DAG are strongly connected components SCC d ijk . Edges of this DAG (purple arrows in Fig. 1) denote solving dependencies between the SCC d ijk and may cross subparts G d ij boundaries. In the following, linkedSCC(G d ij ) refers to the operation that obtains (the subpart of) this DAG from the dth D-M decomposition of CC d ij that corresponds to the given subpart. The third loop consists in trying to iteratively (in the DAGdependencies induced order) find numerical over-constraints in each SCC d ijk or, when it is solvable, propagate its solution to other blocks. Since blocks are strongly connected, there is only one potential solution to each block unless it contains only variables, and this latter case can only be encountered in a third subpart G d i3 . The process works only the top-level of the DAG (red blocks in the figure) because these blocks equations do not use variables from other blocks.

Algorithm 2 1 : 4 :

 214 findRC: Numerical analysis of G d ij subpart of CC d i Require: CC d i and G d ij Ensure: Boolean continue and updated CC d i [SCC d1 ij1 , • • • , SCC d1 ijN ] ←linkedSCC(while [SCC dm ij1 , • • • , SCC dm ijN ] ̸ = ∅ do

Fig. 2 .

 2 Fig. 2. Block analysis of linear systems: (a) overall detection process, (b) QR factorization with column pivoting.

Fig. 3 .

 3 Fig. 3. Block analysis of non-linear systems. Phase I: Over-constraints detection.

Fig. 4 .

 4 Fig. 4. Phase II: Distinguishing redundant and conflicting constraints.

Fig. 7 .

 7 Fig. 7. Results of the sketching after removing conflicting constraints with N rows = 4: (a) initial glass, (b 1 ) configuration 1 and minimization of the shape variation, (b 2 ) configuration 1 and minimization of the area of the final surface, (c 1 ) configuration 3 and minimization of the shape variation, (b 2 ) configuration 3 and minimization of the area of the final surface.

  DoFs and if all the subsystems after decomposition satisfy DoCs ≤ DoFs.

	Definition 4. A geometric constraints system G is structurally
	over-constrained if there exists a subsystem satisfying DoCs >
	DoFs. Structural over-constraints are constraints that transform
	a structurally over-constrained system into a structurally well-
	constrained system when they are removed.
	Definition 2. The degree of freedom DoC (e) of a geometric con-
	straint e is the number of independent equations needed to repre-
	sent it. For instance, distance constraints have one DoC in 2D and
	3D. For a geometric constraints system G with a set E of constraints, the degree of freedom of all the constraints is DoCs = ∑

e∈E DoC (e). Definition 3. A geometric constraints system G is structurally wellconstrained if G satisfies DoCs =

either redundant or conflicting with respect to G b , and if card(E b ) ≥ card(E o ), then E b is a set of basis constraints and E o is a set of numerical over-constraints.

Table 1

 1 Comparison between our algorithm and Moinet's approach on the Double-Banana testing case.

	Method	Witness	Over-constraint	df	cond
	Our	Initial sketch	e9	0.53/45	16.97
	Moinet et al.	Initial sketch	e18	4.38/32	73.33

Table 2

 2 Typology of constraints and equations involved in the description of the 3D glass sketching example.

	Constraint	Equations	Type	Component	Constraint	Equations	Type	Component
	4	1-3	Linear	1	12	30-32	Non-linear	2
	2	4-6	Linear	2	13	33-35	Linear	2
	1	7-9	Linear	1	14	36-38	Non-linear	2
	3	10-12	Linear	2	15	39-41	Linear	1
	5	13	Non-linear	1	16	42-44	Non-linear	1
	6	14	Non-linear	2	17	45-47	Linear	1
	7	15-17	Linear	1	18	48-50	Non-linear	1
	8	18-20	Non-linear	1	19	51-53	Linear	2
	9	21-23	Linear	1	20	54-56	Non-linear	2
	10	24-26	Non-linear	1	21	57-59	Linear	2
	11	27-29	Linear	2	22	60-62	Non-linear	2

Table 3

 3 Status of the distance and position constraints (0 to remove and 1 to keep) to solve the 9 over-constrained configurations.

	Config.	DIS(P1,P4)	DIS(P2,P3)	POS(P1)	POS(P2)	POS(P3)	POS(P4)

Table 4

 4 Evaluation of the 9 configurations with N rows = 4.

	Config.	Minimization of G 1 (X )		Minimization of G 2 (X )
		dg 1	df	cond	dg 2	df	cond
	1	13801.04	0.10000	2.8654e19	96733.72	0.10000	1.4272e18
	2	17990.88	0.10182	8.3172e18	95225.05	0.28157	4.3894e17
	3	15459.52	0.10179	1.5071e19	94483.08	0.22880	4.9533e17
	4	12265.51	0.10975	8.8857e18	89924.13	0.22806	3.9399e18
	5	10970.98	0.10971	3.9852e19	86879.47	0.25225	8.2501e18
	6	15826.68	0.12607	3.4260e18	76878.26	0.68278	1.6567e18
	7	12936.45	0.10465	3.8205e18	76167.99	0.62820	3.9842e17
	8	13889.18	0.11385	2.5681e18	78657.81	0.59160	4.1485e18
	9	14883.21	0.11720	1.2523e18	74351.81	0.71765	6.7658e16

Table 5

 5 Evaluation of the 9 configurations with N rows = 5.

	Config.	Minimization of G 1 (X )		Minimization of G 2 (X )
		dg 1	df	cond	dg 2	df	cond
	1	11266.93	0.10000	4.0149e17	85121.36	0.10000	3.3441e17
	2	14719.05	0.10280	4.6031e17	86295.47	0.25034	6.5355e19
	3	12506.55	0.10277	1.7748e19	85190.96	0.20076	1.0972e18
	4	9944.87	0.11452	1.7903e18	79428.31	0.20592	1.7041e18
	5	8799.29	0.11448	6.1454e17	77800.57	0.22919	1.0218e18
	6	12561.66	0.13935	4.1681e18	69603.16	0.76646	8.5100e16
	7	10441.11	0.10684	1.0862e18	69502.72	0.70009	2.3460e18
	8	11134.09	0.12097	2.5394e18	71465.72	0.65901	1.5773e18
	9	11877.59	0.12601	1.3790e19	67661.55	0.80372	8.4472e17

Acknowledgment

The authors are grateful to the China Scholarship Council (No. 201406090176) for supporting this research. Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org) Powered by TCPDF (www.tcpdf.org)