SINTEF / Splipy

Spline modelling made easy.
GNU General Public License v3.0
102 stars 18 forks source link

Raise Order should not use interpolation #14

Closed VikingScientist closed 4 years ago

VikingScientist commented 7 years ago

Raise order is currently building the correct knot vector followed by Greville-point interpolation. For large models (with small element sizes) this ends up being numerically unstable and numpy is complaining about singular matrices. To trigger this, consider the cylinder.py method in meshscripts and call this with the following arguments

python3 cylinder.py --nel-height 3 --nel-circ 5 --nel-bndl 5 --height 10

Output:

  File "cylinder.py", line 290, in cylinder
    patch.raise_order(0, 0, 2)
  File "/home/kjetijo/Documents/Projects/geomodeler/splipy/SplineObject.py", line 357, in raise_order
    result = np.tensordot(np.linalg.inv(n), result, axes=(1, self.pardim-1))
  File "/usr/local/lib/python3.5/dist-packages/numpy/linalg/linalg.py", line 526, in inv
    ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
  File "/usr/local/lib/python3.5/dist-packages/numpy/linalg/linalg.py", line 90, in _raise_linalgerror_singular
    raise LinAlgError("Singular matrix")
numpy.linalg.linalg.LinAlgError: Singular matrix

The solution is to use a more explicit formulation for the raise-order which does not require the solution of linear system of equations. This will in general improve the entire library as it will be faster execution time with no drawbacks.

VikingScientist commented 4 years ago

The formulas for this can be found in chapter 4.1.2 of Isogeometric analysis of acoustic scattering by J.V.Venås

VikingScientist commented 4 years ago

Solved in #108