add MLine::discretize, improve decasteljau alogrithm

Geo/CMakeLists.txt

@@ -11,7 +11,6 @@ set(SRC
     GEdgeLoop.cpp GEdgeCompound.cpp GFaceCompound.cpp
   GRegionCompound.cpp GRbf.cpp
     gmshVertex.cpp gmshEdge.cpp gmshFace.cpp gmshRegion.cpp
-    gmshEdgeDiscretize.cpp
    gmshSurface.cpp
     OCCVertex.cpp OCCEdge.cpp OCCFace.cpp OCCRegion.cpp
     discreteEdge.cpp discreteFace.cpp discreteRegion.cpp

Geo/MLine.cpp

@@ -10,6 +10,7 @@
 #include "GaussLegendre1D.h"
 #include "Context.h"
 #include "qualityMeasures.h"
+#include "decasteljau.h"
 
 const JacobianBasis* MLine::getJacobianFuncSpace(int order) const
 {
@@ -79,3 +80,62 @@ void MLineN::getEdgeRep(bool curved, int num, double *x, double *y, double *z, S
   }
   else MLine::getEdgeRep(false, num, x, y, z, n);
 }
+
+void MLine::discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts)
+{
+  ts.clear();
+  ts.push_back(-1);
+  ts.push_back(1);
+  dpts.clear();
+  dpts.push_back(getVertex(0)->point());
+  dpts.push_back(getVertex(1)->point());
+}
+
+void MLine3::discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts)
+{
+  SPoint3 p0 = getVertex(0)->point();
+  SPoint3 p2 = getVertex(1)->point();
+  SPoint3 p1 = getVertex(2)->point() * 2 - (p0  + p2) * 0.5;
+  decasteljau(tol, p0, p1, p2, dpts, ts);
+  for (size_t i = 0; i < ts.size(); ++i)
+    ts[i] = -1 + 2 * ts[i];
+}
+
+void MLineN::discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts)
+{
+  int order = getPolynomialOrder();
+  if (order == 3) {
+    SPoint3 p0 = getVertex(0)->point();
+    SPoint3 p3 = getVertex(1)->point();
+    SPoint3 p1 = p0 * (-5./6) + p3 * (1./3) + getVertex(2)->point() * 3. - getVertex(3)->point() * 1.5;
+    SPoint3 p2 = p0 * (1./3) + p3 * (-5./6) - getVertex(2)->point() * 1.5 + getVertex(3)->point() * 3.;
+    decasteljau(tol, p0, p1, p2, p3, dpts, ts);
+    for (size_t i = 0; i < ts.size(); ++i)
+      ts[i] = -1 + 2 * ts[i];
+    return;
+  }
+  fullMatrix<double> lagNodes(order + 1, 3), bezNodes( order + 1, 3);
+  for (int i = 0; i < order + 1; ++i) {
+    MVertex *v = getVertex(i);
+    lagNodes(i, 0) = v->x();
+    lagNodes(i, 1) = v->y();
+    lagNodes(i, 2) = v->z();
+  }
+  const bezierBasis *bez = BasisFactory::getBezierBasis(TYPE_LIN, order);
+  bez->matrixLag2Bez.mult(lagNodes, bezNodes);
+  std::vector<SPoint3> pts(bezNodes.size1());
+  pts[0][0] = bezNodes(0, 0);
+  pts[0][1] = bezNodes(0, 1);
+  pts[0][2] = bezNodes(0, 2);
+  pts[order][0] = bezNodes(1, 0);
+  pts[order][1] = bezNodes(1, 1);
+  pts[order][2] = bezNodes(1, 2);
+  for (int i = 0; i < order - 1; ++i) {
+    pts[i + 1][0] = bezNodes(i + 2, 0);
+    pts[i + 1][1] = bezNodes(i + 2, 1);
+    pts[i + 1][2] = bezNodes(i + 2, 2);
+  }
+  decasteljau(tol, pts, dpts, ts);
+    for (size_t i = 0; i < ts.size(); ++i)
+      ts[i] = -1 + 2 * ts[i];
+}

Geo/MLine.h

@@ -96,6 +96,7 @@ class MLine : public MElement {
     return SPoint3(0, 0, 0);
   }
   virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
+  virtual void discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts);
 };
 
 /*
@@ -149,6 +150,7 @@ class MLine3 : public MLine {
   {
     num < 2 ? MLine::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
+  virtual void discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts);
 };
 
 /*
@@ -208,6 +210,7 @@ class MLineN : public MLine {
   {
     num < 2 ? MLine::getNode(num, u, v, w) : MElement::getNode(num, u, v, w);
   }
+  virtual void discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts);
 };
 
 struct compareMLinePtr {

Geo/gmshEdge.cpp

@@ -12,6 +12,7 @@
 #include "GeoInterpolation.h"
 #include "GmshMessage.h"
 #include "Context.h"
+#include "decasteljau.h"
 
 gmshEdge::gmshEdge(GModel *m, Curve *edge, GVertex *v1, GVertex *v2)
   : GEdge(m, edge->Num, v1, v2), c(edge)
@@ -406,3 +407,127 @@ void gmshEdge::writeGEO(FILE *fp)
   if(meshAttributes.reverseMesh)
     fprintf(fp, "Reverse Line {%d};\n", tag());
 }
+
+static inline SPoint3 curveGetPoint(Curve *c, int i)
+{
+  Vertex *v;
+  List_Read(c->Control_Points, i , &v);
+  return SPoint3(v->Pos.X, v->Pos.Y, v->Pos.Z);
+}
+
+void gmshEdge::discretize(double tol, std::vector<SPoint3> &pts, std::vector<double> &ts)
+{
+  switch(c->Typ) {
+    case MSH_SEGM_LINE :
+      {
+        int NPt = List_Nbr(c->Control_Points);
+        pts.resize(NPt);
+        ts.resize(NPt);
+        for (int i = 0; i < NPt; ++i) {
+          pts[i]= curveGetPoint(c, i);
+          ts[i] = i / (double) (NPt - 1);
+        }
+        return;
+      }
+    case MSH_SEGM_BEZIER :
+      {
+        int NbCurves = (List_Nbr(c->Control_Points) - 1) / 3;
+        for (int iCurve = 0; iCurve < NbCurves; ++iCurve) {
+          double t1 = (iCurve) / (double)(NbCurves);
+          double t2 = (iCurve+1) / (double)(NbCurves);
+          SPoint3 pt[4];
+          for(int i = 0; i < 4; i++) {
+            pt[i] = curveGetPoint(c, iCurve * 3 + i);
+          }
+          std::vector<double> lts;
+          std::vector<SPoint3> lpts;
+          decasteljau(tol, pt[0], pt[1], pt[2], pt[3], lpts, lts);
+          for (size_t i = (iCurve == 0 ? 0 : 1); i < lpts.size(); ++i) {
+            pts.push_back(lpts[i]);
+            ts.push_back(t1 + lts[i] * (t2 - t1));
+          }
+        }
+        break;
+      }
+    case MSH_SEGM_BSPLN:
+      {
+        bool periodic = (c->end == c->beg);
+        int NbControlPoints = List_Nbr(c->Control_Points);
+        int NbCurves = NbControlPoints + (periodic ? -1 : 1);
+        SPoint3 pt[4];
+        for (int iCurve = 0; iCurve < NbCurves; ++iCurve) {
+          double t1 = (iCurve) / (double)(NbCurves);
+          double t2 = (iCurve+1) / (double)(NbCurves);
+          for(int i = 0; i < 4; i++) {
+            int k;
+            if (periodic) {
+              k = (iCurve - 1 + i) % (NbControlPoints - 1);
+              if (k < 0)
+                k += NbControlPoints - 1;
+            }
+            else {
+              k = std::max(0, std::min(iCurve - 2 + i, NbControlPoints -1));
+            }
+            pt[i] = curveGetPoint(c, k);
+          }
+          SPoint3 bpt[4] = {
+            (pt[0] + pt[1] * 4 + pt[2]) * (1./6.),
+            (pt[1] * 2 + pt[2]) * (1./3.),
+            (pt[1] + pt[2] * 2) * (1./3.),
+            (pt[1] + pt[2] * 4 + pt[3]) * (1./6.)
+          };
+          std::vector<double> lts;
+          std::vector<SPoint3> lpts;
+          decasteljau(tol, bpt[0], bpt[1], bpt[2], bpt[3], lpts, lts);
+          for (size_t i = (iCurve == 0 ? 0 : 1); i < lpts.size(); ++i) {
+            pts.push_back(lpts[i]);
+            ts.push_back(t1 + lts[i] * (t2 - t1));
+          }
+        }
+        break;
+      }
+    case MSH_SEGM_SPLN:
+      {
+        int NbCurves = List_Nbr(c->Control_Points) - 1;
+        SPoint3 pt[4];
+        for (int iCurve = 0; iCurve < NbCurves; ++iCurve) {
+          double t1 = (iCurve) / (double)(NbCurves);
+          double t2 = (iCurve+1) / (double)(NbCurves);
+          pt[1] = curveGetPoint(c, iCurve);
+          pt[2] = curveGetPoint(c, iCurve + 1);
+          if(iCurve == 0) {
+            if(c->beg == c->end)
+              pt[0] = curveGetPoint(c, NbCurves - 1);
+            else
+              pt[0] = SPoint3(pt[1] * 2 - pt[2]);
+          }
+          else
+            pt[0] = curveGetPoint(c, iCurve - 1);
+          if(iCurve == NbCurves - 1) {
+            if(c->beg == c->end)
+              pt[3] = curveGetPoint(c, 1);
+            else
+              pt[3] = SPoint3(2 * pt[2] - pt[1]);
+          }
+          else
+            pt[3] = curveGetPoint(c, iCurve + 2);
+          SPoint3 bpt[4] = {
+            pt[1],
+            (pt[1] * 6 + pt[2] - pt[0]) * (1./6.),
+            (pt[2] * 6 - pt[3] + pt[1]) * (1./6.),
+            pt[2]
+          };
+          std::vector<double> lts;
+          std::vector<SPoint3> lpts;
+          decasteljau(tol, bpt[0], bpt[1], bpt[2], bpt[3], lpts, lts);
+          for (size_t i = (iCurve == 0 ? 0 : 1); i < lpts.size(); ++i) {
+            pts.push_back(lpts[i]);
+            ts.push_back(t1 + lts[i] * (t2 - t1));
+          }
+        }
+        break;
+      }
+    default :
+      Msg::Fatal("not implemented");
+  }
+}

Geo/gmshEdgeDiscretize.cpp

-#include <cstdio>
-#include <cmath>
-#include <vector>
-#include "SPoint3.h"
-#include "SVector3.h"
-#include "GEdge.h"
-#include "gmshEdge.h"
-#include "Geo.h"
-
-class discreteList {
-  std::vector<std::pair<SPoint3, double> > _pts;
-  std::vector<int> _next;
-  public:
-  int insertPoint(int pos, const SPoint3 &pt, double t) {
-    _pts.push_back(std::make_pair(pt, t));
-    _next.push_back(_next[pos + 1]);
-    _next[pos + 1] = _pts.size() - 1;
-    return _pts.size() - 1;
-  }
-  void sort(std::vector<SPoint3> &spts, std::vector<double> &ts) {
-    spts.clear();
-    spts.reserve(_pts.size());
-    ts.clear();
-    ts.reserve(_pts.size());
-    for (int p = _next[0]; p != -1; p = _next[p + 1]) {
-      spts.push_back(_pts[p].first);
-      ts.push_back(_pts[p].second);
-    }
-  }
-  discreteList() {
-    _next.push_back(-1);
-  }
-};
-
-static void decasteljau(double tol, discreteList &discrete, int pos, const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2, const SPoint3 &p3, double t0, double t3)
-{
-  SVector3 d30 = p3 - p0;
-  SVector3 d13 = p1 - p3;
-  SVector3 d23 = p2 - p3;
-  SVector3 d130 = crossprod(d13, d30);
-  SVector3 d230 = crossprod(d23, d30);
-  double d = std::max(dot(d130, d130), dot(d230, d230));
-  double l2 = dot(d30, d30);
-
-  if(d < tol * tol * l2) {
-    return;
-  }
-
-  SPoint3 p01((p0 + p1) * 0.5);
-  SPoint3 p12((p1 + p2) * 0.5);
-  SPoint3 p23((p2 + p3) * 0.5);
-  SPoint3 p012((p01 + p12) * 0.5);
-  SPoint3 p123((p12 + p23) * 0.5);
-  SPoint3 p0123((p012 + p123) * 0.5);
-  double t0123 = 0.5 * (t0 + t3);
-  int newpos = discrete.insertPoint(pos, p0123, t0123);
-
-  decasteljau(tol, discrete, pos, p0, p01, p012, p0123, t0, t0123);
-  decasteljau(tol, discrete, newpos, p0123, p123, p23, p3, t0123, t3);
-}
-
-static int discretizeBezier(double tol, discreteList &discrete, int pos, const SPoint3 pt[4], double t0, double t3, bool insertFirstPoint)
-{
-  if (insertFirstPoint)
-    pos = discrete.insertPoint(pos, pt[0], t0);
-  int newp = discrete.insertPoint(pos, pt[3], t3);
-  decasteljau(tol, discrete, pos, pt[0], pt[1], pt[2], pt[3], t0, t3);
-  return newp;
-}
-
-static int discretizeBSpline(double tol, discreteList &discrete, int pos, const SPoint3 pt[4], double t0, double t3, bool insertFirstPoint)
-{
-  SPoint3 bpt[4] = {
-    SPoint3((pt[0] + 4 * pt[1]  + pt[2]) * (1./6.)),
-    SPoint3((2 * pt[1] + pt[2]) * (1./3.)),
-    SPoint3((pt[1] + 2 * pt[2]) * (1./3.)),
-    SPoint3((pt[1] + 4 * pt[2] + pt[3]) * (1./6.))
-  };
-  return discretizeBezier(tol, discrete, pos, bpt, t0, t3, insertFirstPoint);
-}
-
-static int discretizeCatmullRom(double tol, discreteList &discrete, int pos, const SPoint3 pt[4], double t0, double t3, bool insertFirstPoint)
-{
-  SPoint3 bpt[4] = {
-    pt[1],
-    SPoint3(( 6 * pt[1] + pt[2] - pt[0]) * (1./6.)),
-    SPoint3(( 6 * pt[2] - pt[3] + pt[1]) * (1./6.)),
-    pt[2]
-  };
-  return discretizeBezier(tol, discrete, pos, bpt, t0, t3, insertFirstPoint);
-}
-
-static inline SPoint3 curveGetPoint(Curve *c, int i)
-{
-  Vertex *v;
-  List_Read(c->Control_Points, i , &v);
-  return SPoint3(v->Pos.X, v->Pos.Y, v->Pos.Z);
-}
-
-static void discretizeCurve(Curve *c, double tol, std::vector<SPoint3> &pts, std::vector<double> &ts)
-{
-  discreteList discrete;
-  switch(c->Typ) {
-    case MSH_SEGM_LINE :
-      {
-        int NPt = List_Nbr(c->Control_Points);
-        pts.resize(NPt);
-        ts.resize(NPt);
-        for (int i = 0; i < NPt; ++i) {
-          pts[i]= curveGetPoint(c, i);
-          ts[i] = i / (double) (NPt - 1);
-        }
-        return;
-      }
-    case MSH_SEGM_BEZIER :
-      {
-        int back = -1;
-        int NbCurves = (List_Nbr(c->Control_Points) - 1) / 3;
-        for (int iCurve = 0; iCurve < NbCurves; ++iCurve) {
-          double t1 = (iCurve) / (double)(NbCurves);
-          double t2 = (iCurve+1) / (double)(NbCurves);
-          SPoint3 pt[4];
-          for(int i = 0; i < 4; i++) {
-            pt[i] = curveGetPoint(c, iCurve * 3 + i);
-          }
-          back = discretizeBezier(tol, discrete, back, pt, t1, t2, iCurve == 0);
-        }
-        break;
-      }
-    case MSH_SEGM_BSPLN:
-      {
-        int back = -1;
-        bool periodic = (c->end == c->beg);
-        int NbControlPoints = List_Nbr(c->Control_Points);
-        int NbCurves = NbControlPoints + (periodic ? -1 : 1);
-        SPoint3 pt[4];
-        for (int iCurve = 0; iCurve < NbCurves; ++iCurve) {
-          double t1 = (iCurve) / (double)(NbCurves);
-          double t2 = (iCurve+1) / (double)(NbCurves);
-          for(int i = 0; i < 4; i++) {
-            int k;
-            if (periodic) {
-              k = (iCurve - 1 + i) % (NbControlPoints - 1);
-              if (k < 0)
-                k += NbControlPoints - 1;
-            }
-            else {
-              k = std::max(0, std::min(iCurve - 2 + i, NbControlPoints -1));
-            }
-            pt[i] = curveGetPoint(c, k);
-          }
-          back = discretizeBSpline(tol, discrete, back, pt, t1, t2, iCurve == 0);
-        }
-        break;
-      }
-    case MSH_SEGM_SPLN:
-      {
-        int NbCurves = List_Nbr(c->Control_Points) - 1;
-        SPoint3 pt[4];
-        int back = -1;
-        for (int iCurve = 0; iCurve < NbCurves; ++iCurve) {
-          double t1 = (iCurve) / (double)(NbCurves);
-          double t2 = (iCurve+1) / (double)(NbCurves);
-          pt[1] = curveGetPoint(c, iCurve);
-          pt[2] = curveGetPoint(c, iCurve + 1);
-          if(iCurve == 0) {
-            if(c->beg == c->end)
-              pt[0] = curveGetPoint(c, NbCurves - 1);
-            else
-              pt[0] = SPoint3(pt[1] * 2 - pt[2]);
-          }
-          else
-            pt[0] = curveGetPoint(c, iCurve - 1);
-          if(iCurve == NbCurves - 1) {
-            if(c->beg == c->end)
-              pt[3] = curveGetPoint(c, 1);
-            else
-              pt[3] = SPoint3(2 * pt[2] - pt[1]);
-          }
-          else
-            pt[3] = curveGetPoint(c, iCurve + 2);
-          back = discretizeCatmullRom(tol, discrete, back, pt, t1, t2, iCurve == 0);
-        }
-        break;
-      }
-    default :
-      Msg::Fatal("not implemented");
-  }
-  discrete.sort(pts, ts);
-}
-
-void gmshEdge::discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts)
-{
-  discretizeCurve(c, tol, dpts, ts);
-}

Numeric/CMakeLists.txt

@@ -30,6 +30,7 @@ set(SRC
   GaussJacobi1D.cpp
   HilbertCurve.cpp
  robustPredicates.cpp
+ decasteljau.cpp
   mathEvaluator.cpp
   Iso.cpp
 )

Numeric/decasteljau.cpp

+#include "decasteljau.h"
+#include "SPoint3.h"
+#include "SVector3.h"
+
+class discreteList {
+  std::vector<std::pair<SPoint3, double> > _pts;
+  std::vector<int> _next;
+  public:
+  int insertPoint(int pos, const SPoint3 &pt, double t) {
+    _pts.push_back(std::make_pair(pt, t));
+    _next.push_back(_next[pos + 1]);
+    _next[pos + 1] = _pts.size() - 1;
+    return _pts.size() - 1;
+  }
+  void sort(std::vector<SPoint3> &spts, std::vector<double> &ts) {
+    spts.clear();
+    spts.reserve(_pts.size());
+    ts.clear();
+    ts.reserve(_pts.size());
+    for (int p = _next[0]; p != -1; p = _next[p + 1]) {
+      spts.push_back(_pts[p].first);
+      ts.push_back(_pts[p].second);
+    }
+  }
+  discreteList() {
+    _next.push_back(-1);
+  }
+};
+
+static double sqDistPointSegment(const SPoint3 &p, const SPoint3 &s0, const SPoint3 &s1)
+{
+  SVector3 d(s1 - s0);
+  SVector3 d0(p - s0);
+  SVector3 d1(p - s1);
+  double dn2 = crossprod(d, d0).normSq();
+  double dt2 = std::max(0., std::max(-dot(d, d0), dot(d, d1)));
+  dt2 *= dt2;
+  return (dt2 + dn2) / d.normSq();
+}
+
+static void decasteljau(double tol, discreteList &discrete, int pos, const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2, double t0, double t2)
+{
+  if(sqDistPointSegment(p1, p0, p2) < tol * tol)
+    return;
+  SPoint3 p01((p0 + p1) * 0.5);
+  SPoint3 p12((p1 + p2) * 0.5);
+  SPoint3 p012((p01 + p12) * 0.5);
+  double t012 = 0.5 * (t0 + t2);
+  int newpos = discrete.insertPoint(pos, p012, t012);
+  decasteljau(tol, discrete, pos, p0, p01, p012, t0, t012);
+  decasteljau(tol, discrete, newpos, p012, p12, p2, t012, t2);
+}
+
+void decasteljau(double tol, const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2, std::vector<SPoint3> &pts, std::vector<double> &ts)
+{
+  discreteList discrete;
+  discrete.insertPoint(-1, p0, 0);
+  discrete.insertPoint(0, p2, 1);
+  decasteljau(tol, discrete, 0, p0, p1, p2, 0., 1);
+  discrete.sort(pts, ts);
+}
+
+static void decasteljau(double tol, discreteList &discrete, int pos, const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2, const SPoint3 &p3, double t0, double t3)
+{
+  if (std::max(sqDistPointSegment(p1, p0, p3), sqDistPointSegment(p2, p0, p3)) < tol * tol)
+    return;
+  SPoint3 p01((p0 + p1) * 0.5);
+  SPoint3 p12((p1 + p2) * 0.5);
+  SPoint3 p23((p2 + p3) * 0.5);
+  SPoint3 p012((p01 + p12) * 0.5);
+  SPoint3 p123((p12 + p23) * 0.5);
+  SPoint3 p0123((p012 + p123) * 0.5);
+  double t0123 = 0.5 * (t0 + t3);
+  int newpos = discrete.insertPoint(pos, p0123, t0123);
+  decasteljau(tol, discrete, pos, p0, p01, p012, p0123, t0, t0123);
+  decasteljau(tol, discrete, newpos, p0123, p123, p23, p3, t0123, t3);
+}
+
+void decasteljau(double tol, const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2, const SPoint3 &p3, std::vector<SPoint3> &pts, std::vector<double> &ts)
+{
+  discreteList discrete;
+  discrete.insertPoint(-1, p0, 0);
+  discrete.insertPoint(0, p3, 1);
+  decasteljau(tol, discrete, 0, p0, p1, p2, p3, 0., 1);
+  discrete.sort(pts, ts);
+}
+
+static void decasteljau(double tol, discreteList &discrete, int pos, const std::vector<SPoint3> &pts, double t0, double te)
+{
+  int order = pts.size() - 1;
+  double dmax2 = 0;
+  for (int i = 1; i < order ; ++i)
+    dmax2 = std::max(dmax2, sqDistPointSegment(pts[i], pts[0], pts[order]));
+  if(dmax2 < tol * tol)
+    return;
+  std::vector<SPoint3> sub0(pts.size());
+  std::vector<SPoint3> sub1(pts);
+  for (int l = 0; l < order + 1; ++l) {
+    sub0[l] = sub1[0];
+    for (int i = 0; i < order - l; ++i)
+      sub1[i] = (sub1[i] + sub1[i + 1]) * 0.5;
+  }
+  double tmid = 0.5 * (t0 + te);
+  int newpos = discrete.insertPoint(pos, sub1[0], tmid);
+  decasteljau(tol, discrete, pos, sub0, t0, tmid);
+  decasteljau(tol, discrete, newpos, sub1, tmid, te);
+}
+
+void decasteljau(double tol, const std::vector<SPoint3> &controlPoints, std::vector<SPoint3> &pts, std::vector<double> &ts)
+{
+  discreteList discrete;
+  discrete.insertPoint(-1, controlPoints[0], 0);
+  discrete.insertPoint(0, controlPoints.back(), 1);
+  decasteljau(tol, discrete, 0, controlPoints, 0., 1);
+  discrete.sort(pts, ts);
+}

Numeric/decasteljau.h

+#ifndef _DECASTELJAU_H_
+#define _DECASTELJAU_H_
+#include <vector>
+class SPoint3;
+void decasteljau(double tol, const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2, std::vector<SPoint3> &pts, std::vector<double> &ts);
+void decasteljau(double tol, const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2, const SPoint3 &p3, std::vector<SPoint3> &pts, std::vector<double> &ts);
+void decasteljau(double tol, const std::vector<SPoint3> &controlPoints, std::vector<SPoint3> &pts, std::vector<double> &ts);
+double sqDistPointSegment(const SPoint3 &p, const SPoint3 &s0, const SPoint3 &s1);
+#endif

wrappers/gmshpy/gmshGeo.i

@@ -83,6 +83,7 @@ namespace std {
 %include "discreteEdge.h"
 %include "discreteVertex.h"
 %include "discreteRegion.h"
+%include "SPoint3.h"
 %include "MElement.h"
 %include "MVertex.h"
 %include "MTriangle.h"
@@ -94,7 +95,6 @@ namespace std {
 %include "MFace.h"
 %include "MPoint.h"
 %include "SVector3.h"
-%include "SPoint3.h"
 %include "SPoint2.h"
 %include "SBoundingBox3d.h"
 %include "Curvature.h"
@@ -141,3 +141,10 @@ namespace std {
 
 }
 
+%extend MElement {
+  SPoint3 pnt(double xi0, double xi1, double xi2) {
+    SPoint3 p;
+    $self->pnt(xi0, xi1, xi2, p);
+    return p;
+  }
+}