If you write all the flexagon invariants as a system of equations, small deviations around legal states can be written as a linear system. For example, the stiffness of a sheet of paper between (x1,y1) and (x2,y2) provides
(x1 - x2)**2 + (y1 - y2)**2 - L**2 == 0
It can be mitigated to
chi2 = (x1 - x2)**2 + (y1 - y2)**2 - L**2 + other constraints...
The derivatives of chi2 with respect to x1 , x2 , y1 , y2 give linear equations. The system of linear equations is a matrix, and eigenvalue / eigenvector expansion of this matrix gives you linear combinations of parameters x1 , x2 , y1 , y2 that are easy or difficult to bend. Eigenvectors represent a basic set of possible directions, and each corresponding eigenvalue tells you how difficult it is to bend in that direction. Larger eigenvalues are more limited.
The problem with the above is that if there are really valid directions, i.e. chi2 derivative of chi2 with respect to p is 0 (the original restriction is absolutely satisfied), the matrix is singular and cannot be inverted to obtain its own system. If you only want to know what absolutely allowed directions are, you can calculate the zero matrix space instead of your eigensystem. However, I suspect (I never played with phlegan) that the "permitted" directions are associated with a slight bend, in which case chi2 is small but non-zero. Then you will look for small but non-zero eigenvalues. Other degrees of freedom are allowed and uninteresting, such as the translation or rotation of the entire object. To turn it into a pure eigensystem problem (generally empty space), add restrictions to a system with arbitrarily small lambda constants:
chi2 += lambda_x * (x1 + x2)**2/4.0 + lambda_y * (y1 + y2)**2/4.0
You will recognize them in your decision because they will change as you change each lambda. (The example above gives a fine lambda_x for translating to x and lambda_y for translating to y.)
From an implementation point of view, you can use any linear algebra software to calculate solutions and check variations with lambdas. I used Python to prototype such a problem (detection in high-energy physics, in which the constraints are measurements, such as “this detector is 3 cm away from this detector”, and chi2 was obtained from the uncertainties “3 cm + - 0.1 cm” ), and then ported the solution to C ++ (BLAS) for production. There was enough linear algebra in the Numpy library for Python (it has BLAS under the hood), although I also used generalized non-linear minimizers in Scipy to debug the matrix solution. The most difficult task is for the indices to line up correctly, which is necessary when you make it as a matrix, and not when you provide an objective function for a general minimizer (because you use variable names instead). This is more of a Matlab or Mathematica problem, so if you are comfortable with one of them, use it instead. This problem will require a lot of trial and error, so use the most interactive system (with a good REPL interface or worksheet / notebook).
It is also useful to draw a graph of connections (graph graph theory, not graph), on which their limitations can be indicated. For me it was a necessary first step before writing down the equations.
It can also help visualize the system by writing a set of functions that take parameter values ( x1 , etc.) and draw a picture using OpenGL (or other three-dimensional grid rendering). This may show you if any restriction is being violated because the mesh tiles will run against each other. It can also help you determine the degrees of freedom represented by each eigenvector: change the parameters by the linear combination represented by the eigenvector, and you will see if it simply translates / rotates, or if it makes some interesting twist or fold.