diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index eb9fba44b4e5a79d3354e42788c6aaf3def7dd4f..2c56e9b75c9551f9b2e5d9bbc867b66b94357e0b 100644
--- a/Geo/CMakeLists.txt
+++ b/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
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index c0859e53588ff2467dca0f8b6c3902f83b905b9e..1299cb4dc62ef89bb876730e5997d2bb6174f09f 100644
--- a/Geo/MLine.cpp
+++ b/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];
+}
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 98d6d240a7ac0dd22df33a965ff39a4a48f91b14..378a4d62ef3e762ef6ea2d83656ad0eb2883d371 100644
--- a/Geo/MLine.h
+++ b/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 {
diff --git a/Geo/gmshEdge.cpp b/Geo/gmshEdge.cpp
index 8bd39eaf55677799351d061b3696d7b03385e5c7..8d4a2fa398e2f3a15daec307842a49262433b21e 100644
--- a/Geo/gmshEdge.cpp
+++ b/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");
+ }
+}
diff --git a/Geo/gmshEdgeDiscretize.cpp b/Geo/gmshEdgeDiscretize.cpp
deleted file mode 100644
index 440ebe4f7c3ef7dc58b9ae85265b048c2a265720..0000000000000000000000000000000000000000
--- a/Geo/gmshEdgeDiscretize.cpp
+++ /dev/null
@@ -1,195 +0,0 @@
-#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);
-}
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index ae0f52d7217b3b82eaedfd4c9d920379bf992e94..8abf9bc3f3749692942d81cd77b6232bae92901b 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -30,6 +30,7 @@ set(SRC
GaussJacobi1D.cpp
HilbertCurve.cpp
robustPredicates.cpp
+ decasteljau.cpp
mathEvaluator.cpp
Iso.cpp
)
diff --git a/Numeric/decasteljau.cpp b/Numeric/decasteljau.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..15e9c0120ccc533e9c2174212340787886561a46
--- /dev/null
+++ b/Numeric/decasteljau.cpp
@@ -0,0 +1,116 @@
+#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);
+}
diff --git a/Numeric/decasteljau.h b/Numeric/decasteljau.h
new file mode 100644
index 0000000000000000000000000000000000000000..7daa3251cf3007935319bf7a24532069857e8ef3
--- /dev/null
+++ b/Numeric/decasteljau.h
@@ -0,0 +1,9 @@
+#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
diff --git a/wrappers/gmshpy/gmshGeo.i b/wrappers/gmshpy/gmshGeo.i
index 49d6ad1e9f250088c7445ab95226b27d789266e3..67e3432a0c0db5030052b50d753090f8d7ce47db 100644
--- a/wrappers/gmshpy/gmshGeo.i
+++ b/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;
+ }
+}