volume_factory.edge_surfaces does not work?

[UnaiSan]

This function fails to recover the boundaries. The test for it (rebuild a unit cube) works because the geometry is completely linear, but it fails if we do not have a cube:

(this test fails)

import numpy as np

import splipy
import splipy.curve_factory as cf
import splipy.surface_factory as sf
import splipy.volume_factory as vf

def test_edge_surfaces_six_sides():
        # create the unit cube
        vol = splipy.Volume()
        # modify slightly one edge: umin(w=0) is now a parabola instead of a line
        vol.raise_order(0,1,0)
        vol.controlpoints[-1, 1, 0, 0] += 1

        faces = vol.faces()

        # edge_surface should give back the same modified unit cube
        vol2 = vf.edge_surfaces(faces)

        # check discretization
        assert vol2.order(0) == 2
        assert vol2.order(1) == 3
        assert vol2.order(2) == 2

        assert len(vol2.knots(0)) == 2 # [0, 1]
        assert len(vol2.knots(1)) == 2
        assert len(vol2.knots(2)) == 2

        # check a 5x5x5 evaluation grid
        u = np.linspace(0,1,5)
        v = np.linspace(0,1,5)
        w = np.linspace(0,1,5)
        pt  = vol( u,v,w)
        pt2 = vol2(u,v,w)

        assert np.allclose(pt, pt2) # failing here

The solid volume is the original one and the transparent one the one built from faces.

[UnaiSan]

This modification at least makes that test pass.

In the Coons (bilinear) surface we do

  x(0, v) (1 - u) + x(1, v) u
+ x(u, 0) (1 - v) + x(u, 1) v
- [ x(0, 0) (1 - u) (1 - v) + x(0, 1) (1 - u) v
  + x(1, 0) u (1 - v)       + x(1, 1) u v]

So in a Coons volume we should have

  x(0, v, w) (1 - u) + x(1, v, w) u
+ x(u, 0, w) (1 - v) + x(u, 1, w) v
+ x(u, v, 0) (1 - w) + x(u, v, 1) w
- [ x(0, 0, w) (1 - u) (1 - v) + x(0, 1, w) (1 - u) v
  + x(1, 0, w) u (1 - v)       + x(1, 1, w) u v]
- [ x(u, 0, 0) (1 - v) (1 - w) + x(u, 0, 1) (1 - v) w
  + x(u, 1, 0) v (1 - w)       + x(u, 1, 1) v w]
- [ x(0, v, 0) (1 - u) (1 - w) + x(0, v, 1) (1 - u) w
  + x(1, v, 0) u (1 - w)       + x(1, v, 1) u w]

But this volume does not fulfill the boundary conditions at the edges. With the addition of the trilinear vol4 we got it.

        result.controlpoints += vol4.controlpoints

        Nu, Nv, Nw, d = result.controlpoints.shape

        controlpoints = np.zeros((2, 2, Nw, d))
        controlpoints[0, 0, :] = vol1.controlpoints[0, 0]
        controlpoints[0, 1, :] = vol1.controlpoints[0, -1]
        controlpoints[1, 0, :] = vol1.controlpoints[-1, 0]
        controlpoints[1, 1, :] = vol1.controlpoints[-1, -1]
        vol_u_edges = splipy.Volume(
            splipy.BSplineBasis(),
            splipy.BSplineBasis(),
            result.bases[2],
            controlpoints=controlpoints,
            raw=True,
            rational=result.rational
        )

        controlpoints = np.zeros((Nu, 2, 2, d))
        controlpoints[:, 0, 0] = vol2.controlpoints[:, 0, 0]
        controlpoints[:, 0, 1] = vol2.controlpoints[:, 0, -1]
        controlpoints[:, 1, 0] = vol2.controlpoints[:, -1, 0]
        controlpoints[:, 1, 1] = vol2.controlpoints[:, -1, -1]
        vol_v_edges = splipy.Volume(
            result.bases[0],
            splipy.BSplineBasis(),
            splipy.BSplineBasis(),
            controlpoints=controlpoints,
            raw=True,
            rational=result.rational
        )

        controlpoints = np.zeros((2, Nv, 2, d))
        controlpoints[0, :, 0] = vol3.controlpoints[0, :, 0]
        controlpoints[0, :, 1] = vol3.controlpoints[0, :, -1]
        controlpoints[1, :, 0] = vol3.controlpoints[-1, :, 0]
        controlpoints[1, :, 1] = vol3.controlpoints[-1, :, -1]
        vol_w_edges = splipy.Volume(
            splipy.BSplineBasis(),
            result.bases[1],
            splipy.BSplineBasis(),
            controlpoints=controlpoints,
            raw=True,
            rational=result.rational
        )

        splipy.Volume.make_splines_identical(result, vol_u_edges)
        splipy.Volume.make_splines_identical(result, vol_v_edges)
        splipy.Volume.make_splines_identical(result, vol_w_edges)

        result.controlpoints -= vol_u_edges.controlpoints
        result.controlpoints -= vol_v_edges.controlpoints
        result.controlpoints -= vol_w_edges.controlpoints

        return result

    else:
        raise ValueError("Requires two or six input surfaces")

In this image the original volume is rendered (in ParaView) with white wireframe edges and the rebuilt one with solid orange color.

[VikingScientist]

Good job both spotting the bug and providing a suggested fix :+1:

[UnaiSan]

Thanks for this nice library!

I have just found this link with a description of that solution: The transfinite interpolations and applications to mesh an object, navigate to A transfinite interpolation on an hexahedron from the data of 6 FACES. It mentions a reference to a paper by Gordon and Hall, but it turns out that in the Wikipedia page for transfinite interpolation there is a reference to other paper with theory behind this, Gordon and Thiel (1982), Transfinite mapping and their application to grid generation, doi:10.1016/0096-3003(82)90191-6 .