From 480a60c95bcdc0deafe6d5e6c93170f36c16e381 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marcin=20=C5=81o=C5=9B?= <marcin.los.91@gmail.com>
Date: Tue, 25 Apr 2017 22:29:35 +0200
Subject: [PATCH] Add support for bases with separators

---
 src/ads/basis_data.cpp      |  9 ++++-----
 src/ads/bspline/bspline.cpp | 17 ++++++++++++++++-
 src/ads/bspline/bspline.hpp | 14 +++++++++++++-
 3 files changed, 33 insertions(+), 7 deletions(-)

diff --git a/src/ads/basis_data.cpp b/src/ads/basis_data.cpp
index 5900438a..c93aa318 100644
--- a/src/ads/basis_data.cpp
+++ b/src/ads/basis_data.cpp
@@ -24,8 +24,8 @@ namespace ads {
 
         for (int e = 0; e < elements; ++ e) {
             x[e] = new double[q];
-            double x1 = this->basis.knot[p + e];
-            double x2 = this->basis.knot[p + e + 1];
+            double x1 = this->basis.points[e];
+            double x2 = this->basis.points[e + 1];
             J[e] = 0.5 * (x2 - x1);
 
             for (int k = 0; k < q; ++ k) {
@@ -40,9 +40,8 @@ namespace ads {
                 for (int d = 0; d <= derivatives; ++ d) {
                     b[e][k][d] = new double[p + 1];
                 }
-                for (int i = 0; i < this->basis.dofs_per_element(); ++ i) {
-                    eval_basis_with_derivatives(e + p, x[e][k], this->basis, b[e][k], derivatives, ctx);
-                }
+                int span = find_span(x[e][k], this->basis);
+                eval_basis_with_derivatives(span, x[e][k], this->basis, b[e][k], derivatives, ctx);
             }
         }
     }
diff --git a/src/ads/bspline/bspline.cpp b/src/ads/bspline/bspline.cpp
index 6461a392..a2e2c366 100644
--- a/src/ads/bspline/bspline.cpp
+++ b/src/ads/bspline/bspline.cpp
@@ -23,6 +23,22 @@ basis create_basis(double a, double b, int p, int elements) {
     return {std::move(knot), p};
 }
 
+basis create_basis_C0(double a, double b, int p, int elements) {
+    int points = elements + 1;
+    int knot_size = p * points + 2; // clamped B-spline with C^0 separators
+    knot_vector knot(knot_size);
+
+    knot[0] = a;
+    knot[knot_size - 1] = b;
+
+    for (int i = 0; i < points; ++i) {
+        for (int j = 0; j < p; ++ j) {
+            knot[1 + i * p + j] = lerp(i, elements, a, b);
+        }
+    }
+    return {std::move(knot), p};
+}
+
 int find_span(double x, const basis& b) {
     int low = b.begin_idx();
     int high = b.end_idx();
@@ -137,4 +153,3 @@ std::vector<int> first_nonzero_dofs(const basis& b) {
 
 }
 }
-
diff --git a/src/ads/bspline/bspline.hpp b/src/ads/bspline/bspline.hpp
index 46baa748..1f9d1e6d 100644
--- a/src/ads/bspline/bspline.hpp
+++ b/src/ads/bspline/bspline.hpp
@@ -12,6 +12,16 @@ using knot_vector = std::vector<double>;
 struct basis {
     knot_vector knot;
     int degree;
+    std::vector<double> points;
+
+    basis(knot_vector kn, int degree): knot{ std::move(kn) }, degree{ degree } {
+        points.push_back(knot[0]);
+        for (auto i = 1u; i < knot.size(); ++ i) {
+            if (knot[i] != knot[i - 1]) {
+                points.push_back(knot[i]);
+            }
+        }
+    }
 
     std::size_t knot_size() const {
         return knot.size();
@@ -26,7 +36,7 @@ struct basis {
     }
 
     int elements() const {
-        return dofs() - degree;
+        return points.size() - 1;
     }
 
     int begin_idx() const {
@@ -161,6 +171,8 @@ struct eval_ders_ctx : public basis_eval_ctx {
  */
 basis create_basis(double a, double b, int degree, int elements);
 
+basis create_basis_C0(double a, double b, int degree, int elements);
+
 int find_span(double x, const basis& b);
 
 void eval_basis(int i, double x, const basis& b, double* out, basis_eval_ctx& ctx);