GeomAPI_PointsToBSplineSurface generate a C^0 surface when periodic=true

[lavender-lee]

Rencently, I have been learning the OCC by applying it to my problem. Considering that I have a serial points in grid form, and I want to generate a closed surface by GeomAPI_PointsToBSplineSurface:

The corresponding code as listed follows:

#include <GeomAPI_PointsToBSplineSurface.hxx>
#include <TColgp_Array2OfPnt.hxx>
#include <Geom_BSplineCurve.hxx>
#include <Geom_BSplineSurface.hxx>
#include <vector>
#include <iostream>
#include <fstream>
#include <sstream>

std::vector<double> read_data_from_file(const std::string& filename, int& n_rows, int& n_cols) {
    std::ifstream fin(filename);
    std::vector<double> points;
    std::string line;
    n_rows = 0;
    n_cols = 0;
    bool done = false;

    while (std::getline(fin, line)) {
        std::stringstream str_stream(line);
        std::string token;
        while (std::getline(str_stream, token, ',')) {
            if (!done) {
                n_cols++;
            }
            points.push_back(std::stod(token));
        }
        // count the number of columns done
        done = true;
        ++n_rows;
    }
    fin.close();
    return points;
}


Handle(Geom_BSplineSurface) build_spline_surface_by_grid(const double* X, const double* Y, const double* Z, const int* dims, bool row_wise, double TOL) {
    int n_sections = dims[0];
    int n_points = dims[1];
    int idx;
    TColgp_Array2OfPnt grid_points(1, n_points, 1, n_sections);
    for (int i = 0; i < n_sections; i++) {
        for (int j = 0; j < n_points; j++) {
            idx = row_wise ? i * n_points + j : j * n_sections + i;
            grid_points.SetValue(j + 1, i + 1, gp_Pnt(X[idx], Y[idx], Z[idx]));
        }
    }

    //
    GeomAPI_PointsToBSplineSurface algo;
    Standard_Boolean is_periodic = Standard_True;
    algo.Init(grid_points, Approx_ChordLength, 3, 3, GeomAbs_C2, TOL, is_periodic);
    return algo.Surface();
}

int main() {
    std::string filename_X = "D:/data/blade/blade-X.csv";
    std::string filename_Y = "D:/data/blade/blade-Y.csv";
    std::string filename_Z = "D:/data/blade/blade-Z.csv";

    int n_rows, n_cols;
    std::vector<double> X = read_data_from_file(filename_X, n_rows, n_cols);
    std::vector<double> Y = read_data_from_file(filename_Y, n_rows, n_cols);
    std::vector<double> Z = read_data_from_file(filename_Z, n_rows, n_cols);

    std::cout << "n_rows: " << n_rows << " n_cols: " << n_cols << std::endl;

    int dims[2] = { n_rows, n_cols };
    double TOL = 0.01;
    bool row_wise = true;

    Handle(Geom_BSplineSurface)& splineSurf = build_spline_surface_by_grid(&X[0], &Y[0], &Z[0], dims, row_wise, TOL);

    // iso curve in V direction (periodic curve)
    Handle(Geom_Curve) isoCurve = splineSurf->VIso(0);

    Handle(Geom_BSplineCurve) splineCurve = Handle(Geom_BSplineCurve)::DownCast(isoCurve);

    double first = splineCurve->FirstParameter();
    double last = splineCurve->LastParameter();
    gp_Pnt d0_0 = splineCurve->Value(first);
    gp_Pnt d0_1 = splineCurve->Value(last);

    printf("==============================\n");
    printf("periodic curve domain: {%f, %f}\n", first, last);
    printf("0-th der at u=0: {%f, %f, %f}\n", d0_0.X(), d0_0.Y(), d0_0.Z());
    printf("0-th der at u=1: {%f, %f, %f}\n", d0_1.X(), d0_1.Y(), d0_1.Z());

    gp_Vec d1_0 = splineCurve->DN(first, 1);
    gp_Vec d1_1 = splineCurve->DN(last, 1);
    printf("1-th der at u=0: {%f, %f, %f}\n", d1_0.X(), d1_0.Y(), d1_0.Z());
    printf("1-th der at u=1: {%f, %f, %f}\n", d1_1.X(), d1_1.Y(), d1_1.Z());

    gp_Vec d2_0 = splineCurve->DN(first, 2);
    gp_Vec d2_1 = splineCurve->DN(last, 2);
    printf("2-th der at u=0: {%f, %f, %f}\n", d2_0.X(), d2_0.Y(), d2_0.Z());
    printf("2-th der at u=1: {%f, %f, %f}\n", d2_1.X(), d2_1.Y(), d2_1.Z());
    printf("==============================\n");

    return 0;
}

However, the logger in console shows that the iso-curve (a closed curve, periodic curve) is C^0 continuity (In fact, the document said it should be C^1).

Another check is the partial derivative of this surface:

 double v0 = 0.1234; // iso-line

 gp_Pnt S_0 = splineSurf->Value(0, v0);
 gp_Pnt S_1 = splineSurf->Value(1, v0);
 printf("S(0,%f): {%f, %f, %f}\n", v0, S_0.X(), S_0.Y(), S_0.Z());
 printf("S(1,%f): {%f, %f, %f}\n", v0, S_1.X(), S_1.Y(), S_1.Z());

 gp_Vec Su_0 = splineSurf->DN(0, v0, 1, 0);
 gp_Vec Su_1 = splineSurf->DN(1, v0, 1, 0);
 printf("Su(0,%f): {%f, %f, %f}\n", v0, Su_0.X(), Su_0.Y(), Su_0.Z());
 printf("Su(1,%f): {%f, %f, %f}\n", v0, Su_1.X(), Su_1.Y(), Su_1.Z());

S(0,0.123400): {-9.182221, -0.396830, 107.954338}
S(1,0.123400): {-9.182221, -0.396830, 107.954338}
Su(0,0.123400): {0.787309, -0.281236, -5.289451}
Su(1,0.123400): {0.785311, -0.280523, -5.276029}

---

How to explain it? Bug?

[lavender-lee]

In addtion, if I set TOL=-0.01, the PointsToBSplineSurface will use interpolation algorithm. And the resultant surface is C^1 continuity.
Please see the following screenshot:

[lavender-lee]

blade.zip

The debug data can be found in attachment