diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index 2c56e9b75c9551f9b2e5d9bbc867b66b94357e0b..eb9fba44b4e5a79d3354e42788c6aaf3def7dd4f 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -11,6 +11,7 @@ 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/GEdge.h b/Geo/GEdge.h
index 0e5b5cfeb41f58f7c6fdaf43e9c994aedc3fa8aa..cf7b3cbfe9680c41b7a767480a07621da561db44 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -209,6 +209,7 @@ class GEdge : public GEntity {
void addLine(MLine *line){ lines.push_back(line); }
bool computeDistanceFromMeshToGeometry (double &d2, double &dmax);
+ virtual void discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts) {Msg::Fatal("not implemented");}
};
#endif
diff --git a/Geo/gmshEdge.h b/Geo/gmshEdge.h
index f54bcfb4f8f110c6edb439493d34d6cdec33f2aa..3c92afd14a1a54a14261ffb5ab218caae7ceb880 100644
--- a/Geo/gmshEdge.h
+++ b/Geo/gmshEdge.h
@@ -29,6 +29,7 @@ class gmshEdge : public GEdge {
virtual void resetMeshAttributes();
virtual SPoint2 reparamOnFace(const GFace *face, double epar, int dir) const;
virtual void writeGEO(FILE *fp);
+ void discretize(double tol, std::vector<SPoint3> &dpts, std::vector<double> &ts);
};
#endif
diff --git a/Geo/gmshEdgeDiscretize.cpp b/Geo/gmshEdgeDiscretize.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..440ebe4f7c3ef7dc58b9ae85265b048c2a265720
--- /dev/null
+++ b/Geo/gmshEdgeDiscretize.cpp
@@ -0,0 +1,195 @@
+#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/wrappers/gmshpy/gmshGeo.i b/wrappers/gmshpy/gmshGeo.i
index da667b39972f6398a4d9cad475cac75eb6d23704..49d6ad1e9f250088c7445ab95226b27d789266e3 100644
--- a/wrappers/gmshpy/gmshGeo.i
+++ b/wrappers/gmshpy/gmshGeo.i
@@ -58,6 +58,7 @@ namespace std {
%template(DoubleVector) std::vector<double, std::allocator<double> >;
%template(DoubleVectorVector) std::vector<std::vector<double, std::allocator<double> > >;
%template(StringVector) std::vector<std::string, std::allocator<std::string> >;
+ %template(SPoint3Vector) std::vector<SPoint3, std::allocator<SPoint3> >;
}
%include "GmshConfig.h"
@@ -72,6 +73,8 @@ namespace std {
%include "GPoint.h"
%include "GEntity.h"
%include "GVertex.h"
+%apply std::vector<double> &OUTPUT{std::vector<double> &ts}
+%apply std::vector<SPoint3> &OUTPUT{std::vector<SPoint3> &dpts}
%include "GEdge.h"
%include "GFace.h"
%include "GRegion.h"