Extend interpolation to All SplineTypes using static templated methods

Wheels are awesome, but in case we haven't prepared enough wheels, it is nice to have an option to build splintery dynamically.

requires TAB to SPACE

Some improvements are still required to enhance the unit tests, these include:

• add comparison across types to further reduce hard comparisons.
• add order reduction that fails

Some functions that still need bindings to be used in splinepy:

• Conversions between NURBS and Bezier/ rational Bezier
• Derivative evaluation of Bezier/ rational Bezier
• Output of basis functions for Bezier/ rational Bezier (might not be important anyways)
• Reduce_order function for Bezier/ rational Bezier

For me, pre-commit doesn't work with mirrors-clang-format version 16.0.2. After I downgrade the version to 15.0.7 the problem solved. Do we really need the latest version?
Other developers had the same problems: https://github.com/pre-commit/mirrors-clang-format/issues/23

Running the build on splinepy prints the following warning (Ubuntu, python 3.9)

<path>/site-packages/setuptools/command/easy_install.py:144: EasyInstallDeprecationWarning: easy_install command is deprecated. Use build and pip and other standards-based tools.
  warnings.warn(
<path>/site-packages/setuptools/command/install.py:34: SetuptoolsDeprecationWarning: setup.py install is deprecated. Use build and pip and other standards-based tools.

Is this an issue with my installation or a common problem?

I've realized, that there's some inconsistency.

core func:

std::shared_ptr<PySpline>
Compose(const std::shared_ptr<PySpline>& inner,
        const std::shared_ptr<PySpline>& outer) { ... }

python func:

composed = splinepy_core.compose(self, inner_function)

For compose, I guess you'd "compose" inner_function into the outer_function. This means, core func has the flipped arg naming -- is it correct @jzwar?

Recently, Github started to offer codespaces with 60-120 core hours a month for free. Would this maybe be an opportunity to provide jupyter notebooks to test virtually prior to downloading or installing splinepy? I think this would make a good marketing strategy 😉

I still have to figure out the details, just wanted to throw it in the ring

In my opinion, a branch that does not pass the tests should not be merged into the main branch. How could this code end up on the main branch? The code must pass all tests before it is merged (pull-request) in the main branch.
Is there a bug in the CI?

Since the compiler got a bit more "pedantic", a lot of warnings are thrown during the build of splinepy on the c++ side.

Most of them are relatively unimportant and can be resolved by the compiler, however, it risks overlooking some relavant issues. As is good practice, we should try to get rid of these warnings asap.

import splinepy as sp
import numpy as np
a = sp.Bezier(degrees=[1,1,1],control_points=sp.utils.data.cartesian_product([np.array([0,1]) for _ in range(3)]))
aa = a.copy()
aa.cps += [1,0,0]
b = sp.Multipatch(splines=[a,aa])
# First
v = b.boundary_multi_patch()
v.splines # Works fine
v = b.boundary_multi_patch()
v.splines # Error

The above example produces an error as the shared pointer (member of py_multi_patch) is not converted into the correct type. Maybe @j042 you have an idea what causes this problem?

When executing npz.py in splinepy/io a ModuleNotFoundError will be raised.

This is because of the following line of code:
import numpy as np

Numpy is using json but because of the fact that json.py is another file in the splinepy/io folder, it will therefore execute this file instead.
This leads to the following error message:

Traceback (most recent call last):
  File "...\splinepy\splinepy\io\npz.py", line 11, in <module>
    import numpy as np
  File "...\AppData\Local\python\mu\mu_venv-38-20221219-011508\lib\site-packages\numpy\__init__.py", line 140, in <module>
    from . import core
  File "...\AppData\Local\python\mu\mu_venv-38-20221219-011508\lib\site-packages\numpy\core\__init__.py", line 9, in <module>
    from numpy.version import version as __version__
  File "...\AppData\Local\python\mu\mu_venv-38-20221219-011508\lib\site-packages\numpy\version.py", line 3, in <module>
    from ._version import get_versions
  File "...\AppData\Local\python\mu\mu_venv-38-20221219-011508\lib\site-packages\numpy\_version.py", line 7, in <module>
    import json
  File "...\splinepy\splinepy\io\json.py", line 10, in <module>
    from splinepy import settings
ModuleNotFoundError: No module named 'splinepy'

import splinepy as sp
import numpy as np
sp.utils.data.cartesian_product([np.array([0,1,2]),np.array([0,2,4])])

Error:

  File "<stdin>", line 1, in <module>
  File "/home/zwar/Git/gustav_dev/splinepy/splinepy/utils/data.py", line 312, in cartesian_product
    dim_reducing_views[i][..., i] = arrays[i][idx[: n_arr - i]]
ValueError: could not broadcast input array from shape (3,) into shape (2,)

There seems to be some functionality in our CI, that is deprecated.
Check the warning, e.g., here

build_and_tests (3.8, ubuntu-20.04)
The `set-output` command is deprecated and will be disabled soon. Please upgrade to using 
Environment Files. For more information see: 
https://github.blog/changelog/2022-10-11-github-actions-deprecating-save-state-and-set-output-commands/

Would some of the CI expert please check this out :D

is quite outdated and needs updates and (at least) a round trip test

Again as in gustaf the
pip install splinepy[all]
does not work...

and even worse vedo is not listed as an (optional) dependency.

Konstantin hat heute versucht zu instalieren.

Problem
The routine splinepy.BSpline.interpolate_surface does not work for an unequal number of control points in the two parametric directions.

Example
Please find an example below which illustrates how an easy surface can be interpolated for an equal number of control points, but not for an unequal number.
The code to produce the example is documented below the images.

Equal number of control points (3,3)

Unequal number of control points (3,4)

Unequal number of control points (4,3)

import numpy as np
import splinepy


def interpolate_function(f, samples_per_direction):
    """Interpolate the function f with a given number of samples 
    in each direction on the interval [-1, 1] x [-1, 1]

    Parameters
    ----------
    f : function
        Function to interpolate
    samples_per_direction : list
        Number of samples per direction

    Returns
    -------
    np.ndarray, splinepy.BSpline
        Target points and interpolating BSpline
    """

    # Create dataset
    x = [np.linspace(-1, 1, samples_per_direction[0]), 
        np.linspace(-1, 1, samples_per_direction[1])]
    xx, yy = np.meshgrid(x[0], x[1])

    target_points = np.vstack(
        (xx.flatten(), yy.flatten(), h(xx, yy).flatten())
    ).T

    interpolated_surface = splinepy.BSpline.interpolate_surface(
        target_points,
        size_u=samples_per_direction[0],
        size_v=samples_per_direction[1],
        degree_u=2,
        degree_v=2,
    )
    return target_points, interpolated_surface
    
def h(x, y):
    return x + y

if __name__ == "__main__":
    try:
        import gustaf as gus

        gustaf_available = True
    except ImportError:
        gustaf_available = False


    for sample_sizes in [[3, 3], [3, 4], [4, 3]]:
        target_points, interpolated_surface = interpolate_function(h, sample_sizes)
        
        if gustaf_available:
            bspline = gus.BSpline(**interpolated_surface.todict())
            gus.show.show_vedo([bspline])

Solution
In cpp/splinepy/fitting/fitting.cpp the points are ordered according to their x-coordinate in line 85. There is an index mistake here, which affects the computation only for size_u != size_v.
I am happy to provide a PR with a solution.

initializing a spline with a spline will takes backend spline, even if the spline has local changes.
Undesired behavior

else breaks all the index assumptions at cpp side.

During clueless and clumsy playing around with splinepy (0.0.26 from pypi), I triggerd several segfaults, when passing invalid arguments to functions.
I found some more, but have not saved the programs that crashed...
I used larry to run those programs, with anaconda 2021.11 (py3.9).

Invalid parametric coordinate when calling jacobian:

import numpy as np
import splinepy

bspline1 = splinepy.BSpline(
        degrees=[3],
        knot_vectors=[[0, 0, 0, 0, 1, 1, 1, 1]],
        control_points=np.array([[0.0, 0.0], [0.2, 3.0], [1.2, 2.0], [2.0, 2.0]]),
)

print(bspline1.jacobian([[1.0]]))  # OK
print(bspline1.sample(10))  # OK
print(bspline1.jacobian([[1.1]]))  # Segfault

again, but with unclamped knot vector, here also sample crashes:

import numpy as np
import splinepy

bspline1 = splinepy.BSpline(
        degrees=[3],
        # Create an unclamped, uniform knot vector
        knot_vectors=[[0, 1, 2, 3, 4, 5, 6, 7]],
        control_points=np.array([[0.0, 0.0], [0.2, 3.0], [1.2, 2.0], [2.0, 2.0]]),
)

print(bspline1.knot_vector_bounds)

print(bspline1.jacobian([[3.5]]))  # OK
print(bspline1.jacobian([[4.0]]))  # segfault
print(bspline1.sample(10))  # Segfault 

There is also a segfault when calling bspline1.show() either in a jupyter notebook or in an environment that does not have a physical GPU but only MESA. Not sure if you can detect that and simply write an error message or if that has to be fixed in vedo?

Would be nice, to have vector graphics exports for publications etc. svg has a simple syntax and can be translated into all kinds of other formats, e.g., eps pdf using Inkscape for example.

I think it has something to do with python __setitem__ syntax, but it'd be nice if we can avoid it

For some applications it would be useful to have a consistent enumeration of boundaries. The numbering system is somehow arbitrary and changes depending on the FE-System in use. As we are theoretically independent of the dimension, I would suggest the following numbering system

for $j\in{0,\dots, n_{dim}-1}$ axis dimensions. This is applicable to higher dimensions as well

This would require to change the layout of the extract_boundary routine

We would like to have a thin logger functionality for splinepy, similar to what has been implemented in gustaf.

The edge cases seem to be redundant and result in tangled splines. I don't know where they came from but we might need to check their consistency!

Concerns all splines 12-end in the create_tile-function

We should extend the Readme.md file and add some information for newcomers.

Ideas can be collected in this thread:

• CI & CD status with more detail
• Gustaf Logo
• Link to documentation
• Link to thirdparty librarires
• Feature list
• Citing information
• Link to contribution info
• Some details

After all, it is the very first thing people will check out

We need some function normalize_weights, that sets the biggest weight to 1. All multiplications of the weight vector with a scalar do not change the representation of the spline geometry.

$$ C(x) = \frac{\sum_i w_i C_i B_i(x)}{\sum_j w_j B_j(x)} = \frac{\sum_i (\alpha w_i) C_i B_i(x)}{\sum_j (\alpha w_j) B_j(x)} $$

However, in order to perform compositions between splines, the weighted control points must fit the ND unit cube.

This spline

rational = sp.RationalBezier(
    degrees=[2],
    control_points=[[0,0],[0.5,1],[1,0]],
    weights=[1,2,1]
)

should work as this on:

rational = sp.RationalBezier(
    degrees=[2],
    control_points=[[0,0],[0.5,1],[1,0]],
    weights=[0.5,1,0.5]
)

however, the former raises an exception...

Plotting a function on a spline is very complicated at the moment, as every time the resolution and on case need to be dealt with. It is not possible to just plot the same function on a Multipatch. That makes it very hard for new users to understand.

I think it would be nice to have a layout like this

def foo(spline, queries):
  # Function in physical space
  x = spline.evaluate(queries)
  return x*x

def bar(spline, queries):
  # Function defined in parametric space
  return queries*queries

Potential Enhancements

Features and good/cool ideas can be documented and discussed here. This main part will be updated regularly.
last update: 06.10.2022

• N-dim bezier operations
• Multithread executions for all queries
• Dynamic cmake build
• Explicit instantiation of some template classes for faster compilation (develop)
• One example for each functionalities.
• static fitting functions (#24)
• Enhance IO, especially mfem, (#31)
• Boundary spline extraction, especially useful for proximity queries.
• cmake extension - real packaging for cpp use

instead of x0, x1, .... It seems to be hardcoded for basic control points of 1-4 Dim splines.

When the axes of two patches are rotated by pi radians, an incorrect orientation is calculated. In the following example, an alignment of [1, 0] and an orientation of [0, 0] are calculated. However, in my opinion, it should result alignment of [1, 0] and an orientation of [1, 1].

The following test case fails because of the wrong orientations:

def test_orientation_rotated_pi_radians(self):
        # Init Splines to be tested
        a_s = c.splinepy.Bezier(
            degrees=[1, 1], control_points=[[1, 1], [1, 2], [2, 1], [2, 2]]
        )
        b_s = c.splinepy.Bezier(
            degrees=[1, 1], control_points=[[2, 3], [1, 3], [2, 2], [1, 2]]
        )


        # Expected Values
        expected_interfaces = [
            [-1, 7, -1, -1],
            [-1, -1, -1, 1],
        ]

        expected_orientations = [
            [0, 1, 1, 3, 1, 0, 1, 1],
        ]

        # Provide connectivity data
        mp = c.splinepy.Multipatch([a_s, b_s])

        # First compare interfaces to what is expected
        interfaces = mp.interfaces
        self.assertTrue(c.np.all(interfaces == expected_interfaces))

        # Then check orientations
        orientations = mp.interface_orientations()
        self.assertTrue(c.np.all(orientations == expected_orientations))

some tiles have Tile at the end and some doesn't. Can they be unified without Tile at the end?

If colors are set differently for the knots and control points via keyword arguments, these will be ignored if the color c is stated. Here is an example:

import splinepy as sp

# Options
options = {
    "control_points": True,
    "control_point_ids": False,
    "control_mesh_c": "k",  # Will be ignored
    "control_point_r": 20,
    "control_point_c": "k",  # Will be ignored
    "c": (0, 240, 156),  # Uses this instead : why?
    "lighting": "off",
}

deformation_function = sp.BSpline(
    degrees=[2, 1],
    knot_vectors=[[0, 0, 0, 0.5, 1, 1, 1], [0, 0, 1, 1]],
    control_points=[
        [0.0, 0.0],
        [0.5, 0.2],
        [1.0, 0.],
        [1.5, 0.2],
        [0.0, 1.0],
        [0.5, 1.2],
        [1.0, 1.0],
        [1.5, 1.2],

    ]
)
deformation_function.show(**options)

Results in:

wrong return type casting in splines/rational_bezier.hpp

let's go with fname, splines

The pipeline docs/minimal_explicit_build_and_docs fails with the error:

Warning, treated as error:
The default value for `navigation_with_keys` will change to `False` in the next release. If you wish to preserve the old behavior for your site, set `navigation_with_keys=True` in the `html_theme_options` dict in your `conf.py` file.Be aware that `navigation_with_keys = True` has negative accessibility implications:https://github.com/pydata/pydata-sphinx-theme/issues/1492

@clemens-fricke : see examples/parallel_evaluation.py

causes error, as it always search for additional knots.

https://tataratat.github.io/splinepy/_generated/splinepy.helpme.integrate.Integrator.html#splinepy.helpme.integrate.Integrator

I tried to export a microstructure to .IGS files format and could not open it in any other software, so I tried to reimport it into splinepy and got the following error:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
[c:\Users\mkofler\OneDrive](file:///C:/Users/mkofler/OneDrive) - TU Wien\Documents\05_Lattice Structure Facade\minimum_working_example.ipynb Cell 1 line 1
     [15](vscode-notebook-cell:/c%3A/Users/mkofler/OneDrive%20-%20TU%20Wien/Documents/05_Lattice%20Structure%20Facade/minimum_working_example.ipynb#W0sZmlsZQ%3D%3D?line=14) microstructure.tiling = [8, 4, 1]
     [17](vscode-notebook-cell:/c%3A/Users/mkofler/OneDrive%20-%20TU%20Wien/Documents/05_Lattice%20Structure%20Facade/minimum_working_example.ipynb#W0sZmlsZQ%3D%3D?line=16) sp.io.iges.export("test.igs",microstructure.create().patches)
---> [18](vscode-notebook-cell:/c%3A/Users/mkofler/OneDrive%20-%20TU%20Wien/Documents/05_Lattice%20Structure%20Facade/minimum_working_example.ipynb#W0sZmlsZQ%3D%3D?line=17) ig = sp.io.iges.load("test.igs")

File [c:\Users\mkofler\miniconda3\envs\splinepy\Lib\site-packages\splinepy\io\iges.py:21](file:///C:/Users/mkofler/miniconda3/envs/splinepy/Lib/site-packages/splinepy/io/iges.py:21), in load(fname)
      8 def load(fname):
      9     """
     10     Read spline in `.iges` form.
     11 
   (...)
     19       Spline Type defined in NAME_TO_TYPE
     20     """
---> 21     return ioutils.dict_to_spline(splinepy_core.read_iges(fname))
ValueError: invalid stoi argument

Minimum working example:

import splinepy as sp
microstructure = sp.microstructure.Microstructure()

base2D = sp.Bezier(
    degrees=[1, 1],
    control_points=[[ 0.,    0.,    0.  ],
                    [ 2.,    0.,    0.  ],
                    [ 0.,    1.5,   0 ],
                    [ 2.,    1.5,   0 ]],
)

microstructure.deformation_function = base2D.create.extruded(extrusion_vector=[0, 0, 1])

microstructure.microtile = getattr(sp.microstructure.tiles, "NutTile3D")().create_tile()[0]
microstructure.tiling = [8, 4, 1]

sp.io.iges.export("test.igs",microstructure.create().patches)
ig = sp.io.iges.load("test.igs")

List of required features for multipatch geometries

• interelement orientation (automatized) for non-structured configurations
• boundary extraction
• field addition (for gismo or other related export types)

Update the import functions with NAME_TO_TYPE constructors, so we can avoid passing dicts:

• IGES
• MFEM
• JSON
• GISMO
• IRIT
• NPZ
• XML CATS-style (is this even still used?)

Minimal example:

import splinepy as sp

# Create a simple spline
a = sp.BSpline([1],[[0,0,1,2,2]],[[0],[1],[2]])
# Remove first knot (optional)
a.remove_knots(0,[1]) # works just fine
# Error 
a.remove_knots(0,[0]) # Segfault
# Similarly
a.remove_knots(0,[2]) # Segfault

In order to describe complex geometries, more than one spline-patch is required. Currently they are stored as a simple list of splines. However, for analysis, more information might be required, such as:

• Connectivity information
• Boundary information
• (optional) material IDs

Having such a class, would facilitate the multipatch output for several IGA-solvers. It could further provide a few functions

• show for multipatch, storing already sampled geometries
• multithread sample-ing
• select_boundary for interactively chosing boundaries based on patch and face IDs

We should also make sure, that intermediate "designs" can be stored in json format, so boundaries can be stored and altered later!

This is related to gismo, mfem and nutils exports and might be of interest to the users of these softwares.

multipatch export for mfem causes std::bad_alloc if you try to read it with mfem. any idea @jzwar ?

Current show_microstructure.py throws an error, because the bezier-extraction is requested for a bezier type and there is no longer a safeguard for this. Error is already fixed on #95

MFEM export for multipatch systems

Using Bezman routins

In a current PR in bezman, MFEM export was restructured to only use corner vertices of the splines, and it provides functionality to identify the connectivity, the edge enumeration with identical knot vectors as well as boundary elements. The use of corner-vertices-only makes it available also for arbitrary spline types (which could be handed on the python side).

Problem with misalignment

Still, the MFEM export is not "safe", as is relies on some strong assumptions. MFEM export requires spline patches to be structured, i.e., neighboring elements need to be ordered in the same manner.
Example:

works:
3 --- 2 3 --- 2
|     | |     |
0-----1 0 --- 1
does not work (knot vectors in eta direction misaligned):
3 --- 2 1 --- 0
|     | |     |
0-----1 2 --- 3

This convention is arbitrary, as knotvectors can still be matching between two faces. Further, parametric axis can be arbitrarily changed within a spline. Starting from a random seed, we can propagate this information and permutate axis, which would allow the mfem export for all "structurable" spline patch systems.

Proposal

I would suggest solving this issue in two steps

• Use bezman's ExtractMFEMInformation routine in splinepy. This could be a fast solution to make splinepy accessible for first tests
• Implement permute functionality in splinepy and solve mfem export in two steps. If the routine succeeds in a first try, 🥳 if not, retrieve the connectivity from MFEM and use the permutation functions starting from a seed to structure mesh, let it fail if mesh non-restructurable... imo, restructuring is not super expensive, and should probably be done in python (changing and - if required - flipping knotvectors, reordering ctps).

The first step could be done very fast, the second one would be great for the next HACK DAYS ™️

@j042 I would love to hear your opinion on this...