Modify Tsai-Wu to work with safety factors

[A-CGray]

Currently, the Tsai-Wu failure evaluation directly returns the value of the quadratic Tsai-Wu criterion:

$$fail = F_1\sigma_{11} + F_2\sigma{22} + F_{11}\sigma_{11}^2 + F_{22}\sigma_{22}^2 + F_{66}\sigma_{12}^2 + F_{12}\sigma_{11}\sigma_{22} = \sigma^T F \sigma + f^T \sigma $$

This has the disadvantage that the failure value does not scale linearly with the stresses and thus the typical approach of applying a safety factor to a failure constraint during optimisation is not valid. In other words, constraining the structure to have a maximum failure value of 1/2 does not mean the structure has a factor of safety of 2.

$$fail(a\sigma) \neq a fail(\sigma)$$

The von-Mises is also a quadratic failure criteria, but we get around the issue above by using a square root:

$$ fail = \frac{\sigma_{vm}}{\sigma_y} = \frac{\sqrt{\sigma_{11}^2 + \sigma_{22}^2 - \sigma_{11}\sigma_{22} + 3\sigma_{12}^2}}{\sigma_y}$$

For Tsai-Wu, we can't simply take a square root because the failure criteria contains quadratic and linear terms and can sometimes be negative. Instead we have to find the factor by which we would need to scale up/down the stresses to be at the failure boundary:

$$fail = \alpha, \quad s. t. (\frac{1}{\alpha}\sigma)^T F (\frac{1}{\alpha}\sigma) + f^T (\frac{1}{\alpha}\sigma) = 1$$

This equation has the roots

$$\alpha = 1/2 \left( f^T \sigma \pm \sqrt{(f^T \sigma)^2 + 4\sigma^T F \sigma}\right)$$

Since the Tsai-Wu criterion is only physical if the F matrix is positive definite, the only positive root is

$$\alpha = 1/2 \left( f^T \sigma + \sqrt{(f^T \sigma)^2 + 4\sigma^T F \sigma}\right)$$

This PR modifies the Tsai-Wu failure criterion calculation in the Orthotropic ply class so that the value returned is the $\alpha$ value from the above equation. Which means that the criterion can safely be used with safety factors and that the OrthotropicPly class returns the von-Mises failure criteria when used with an isotropic material.

I also implemented the necessary derivatives and verified them against FD using testFailureSens:

Testing failure sensitivity for angle = 21.37 
Failure criteria =  1.92497063e+04 
sens[0] FD:  1.53727582e+02 An:  1.53727580e+02 Error:  1.248e-08 
sens[1] FD:  9.22496438e+01 An:  9.22496429e+01 Error:  9.210e-09 
sens[2] FD:  1.12158574e+02 An:  1.12158577e+02 Error: -3.133e-08 
Angle sensitivity FD:  1.08600003e+04 An:  1.08600003e+04 Error:  1.147e-09 
Failure criteria =  1.91585547e+04 
sens[0] FD:  1.54320900e+02 An:  1.54320901e+02 Error: -4.408e-09 
sens[1] FD:  9.00500727e+01 An:  9.00500722e+01 Error:  5.227e-09 
sens[2] FD:  1.12155638e+02 An:  1.12155640e+02 Error: -1.731e-08 
Angle sensitivity FD:  1.11188533e+04 An:  1.11188533e+04 Error:  1.232e-09 
Failure criteria =  1.80378991e+02 
sens[0] FD: -7.78743527e+01 An: -7.78743526e+01 Error:  2.184e-10 
sens[1] FD: -6.49587185e+01 An: -6.49587185e+01 Error:  3.158e-10 
sens[2] FD: -3.75459197e+01 An: -3.75459197e+01 Error:  1.140e-10 
Angle sensitivity FD:  1.29156342e+01 An:  1.29156341e+01 Error:  1.930e-09

This idea is based on https://doc.comsol.com/6.0/doc/com.comsol.help.sme/sme_ug_theory.06.43.html

[A-CGray]

Tests are failing because these changes affect some of the reference values that are tested against, I can update them if the two of you are happy with these changes

[timryanb]

This makes sense to me. Let's update the reference values so the composite tests pass

[timryanb]

Let's give @gjkennedy a chance to review before merging

[gjkennedy]

I do have some concerns with this approach.

1. This is switching from a failure index to a failure load computation. Yes, this now makes the computation scale linearly with load. However, scaling the strain or stress by the load factor would have the same result. How will this change the performance on different optimization problems?
2. Whenever I use the quadratic formula, I like to consider the possibility of subtractive cancellation: https://www.johndcook.com/blog/2018/04/28/quadratic-formula/
3. We should keep the old failure code in place. This is useful for lamina-level computations.

[A-CGray]

@gjkennedy I added back the original Tsai-Wu index as an option but made the strength ratio version the default, is this good with you?

[A-CGray]

The issues that I thought were in these Tsai-Wu sensititivities actually turned out to be in the Von Mises criterion sensitivities (see #212 ) so this PR is good to be merged.