1 // Gmsh - Copyright (C) 1997-2021 C. Geuzaine, J.-F. Remacle
2 //
3 // See the LICENSE.txt file in the Gmsh root directory for license information.
4 // Please report all issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
5
6 #include <cstdio>
7 #include <cmath>
8 #include <vector>
9 #include "SPoint3.h"
10 #include "SVector3.h"
11 #include "GEdge.h"
12 #include "gmshEdge.h"
13 #include "Geo.h"
14
15 class discreteList {
16 std::vector<std::pair<SPoint3, double> > _pts;
17 std::vector<int> _next;
18
19 public:
insertPoint(int pos,const SPoint3 & pt,double t)20 int insertPoint(int pos, const SPoint3 &pt, double t)
21 {
22 _pts.push_back(std::make_pair(pt, t));
23 _next.push_back(_next[pos + 1]);
24 _next[pos + 1] = _pts.size() - 1;
25 return _pts.size() - 1;
26 }
sort(std::vector<SPoint3> & spts,std::vector<double> & ts)27 void sort(std::vector<SPoint3> &spts, std::vector<double> &ts)
28 {
29 spts.clear();
30 spts.reserve(_pts.size());
31 ts.clear();
32 ts.reserve(_pts.size());
33 for(int p = _next[0]; p != -1; p = _next[p + 1]) {
34 spts.push_back(_pts[p].first);
35 ts.push_back(_pts[p].second);
36 }
37 }
discreteList()38 discreteList() { _next.push_back(-1); }
39 };
40
decasteljau(double tol,discreteList & discrete,int pos,const SPoint3 & p0,const SPoint3 & p1,const SPoint3 & p2,const SPoint3 & p3,double t0,double t3)41 static void decasteljau(double tol, discreteList &discrete, int pos,
42 const SPoint3 &p0, const SPoint3 &p1, const SPoint3 &p2,
43 const SPoint3 &p3, double t0, double t3)
44 {
45 SVector3 d30 = p3 - p0;
46 SVector3 d13 = p1 - p3;
47 SVector3 d23 = p2 - p3;
48 SVector3 d130 = crossprod(d13, d30);
49 SVector3 d230 = crossprod(d23, d30);
50 double d = std::max(dot(d130, d130), dot(d230, d230));
51 double l2 = dot(d30, d30);
52
53 if(d < tol * tol * l2) { return; }
54
55 SPoint3 p01((p0 + p1) * 0.5);
56 SPoint3 p12((p1 + p2) * 0.5);
57 SPoint3 p23((p2 + p3) * 0.5);
58 SPoint3 p012((p01 + p12) * 0.5);
59 SPoint3 p123((p12 + p23) * 0.5);
60 SPoint3 p0123((p012 + p123) * 0.5);
61 double t0123 = 0.5 * (t0 + t3);
62 int newpos = discrete.insertPoint(pos, p0123, t0123);
63
64 decasteljau(tol, discrete, pos, p0, p01, p012, p0123, t0, t0123);
65 decasteljau(tol, discrete, newpos, p0123, p123, p23, p3, t0123, t3);
66 }
67
discretizeBezier(double tol,discreteList & discrete,int pos,const SPoint3 pt[4],double t0,double t3,bool insertFirstPoint)68 static int discretizeBezier(double tol, discreteList &discrete, int pos,
69 const SPoint3 pt[4], double t0, double t3,
70 bool insertFirstPoint)
71 {
72 if(insertFirstPoint) pos = discrete.insertPoint(pos, pt[0], t0);
73 int newp = discrete.insertPoint(pos, pt[3], t3);
74 decasteljau(tol, discrete, pos, pt[0], pt[1], pt[2], pt[3], t0, t3);
75 return newp;
76 }
77
discretizeBSpline(double tol,discreteList & discrete,int pos,const SPoint3 pt[4],double t0,double t3,bool insertFirstPoint)78 static int discretizeBSpline(double tol, discreteList &discrete, int pos,
79 const SPoint3 pt[4], double t0, double t3,
80 bool insertFirstPoint)
81 {
82 SPoint3 bpt[4] = {SPoint3((pt[0] + 4 * pt[1] + pt[2]) * (1. / 6.)),
83 SPoint3((2 * pt[1] + pt[2]) * (1. / 3.)),
84 SPoint3((pt[1] + 2 * pt[2]) * (1. / 3.)),
85 SPoint3((pt[1] + 4 * pt[2] + pt[3]) * (1. / 6.))};
86 return discretizeBezier(tol, discrete, pos, bpt, t0, t3, insertFirstPoint);
87 }
88
discretizeCatmullRom(double tol,discreteList & discrete,int pos,const SPoint3 pt[4],double t0,double t3,bool insertFirstPoint)89 static int discretizeCatmullRom(double tol, discreteList &discrete, int pos,
90 const SPoint3 pt[4], double t0, double t3,
91 bool insertFirstPoint)
92 {
93 SPoint3 bpt[4] = {pt[1], SPoint3((6 * pt[1] + pt[2] - pt[0]) * (1. / 6.)),
94 SPoint3((6 * pt[2] - pt[3] + pt[1]) * (1. / 6.)), pt[2]};
95 return discretizeBezier(tol, discrete, pos, bpt, t0, t3, insertFirstPoint);
96 }
97
curveGetPoint(Curve * c,int i)98 static inline SPoint3 curveGetPoint(Curve *c, int i)
99 {
100 Vertex *v;
101 List_Read(c->Control_Points, i, &v);
102 return SPoint3(v->Pos.X, v->Pos.Y, v->Pos.Z);
103 }
104
discretizeCurve(Curve * c,double tol,std::vector<SPoint3> & pts,std::vector<double> & ts)105 static void discretizeCurve(Curve *c, double tol, std::vector<SPoint3> &pts,
106 std::vector<double> &ts)
107 {
108 discreteList discrete;
109 switch(c->Typ) {
110 case MSH_SEGM_LINE: {
111 int NPt = List_Nbr(c->Control_Points);
112 pts.resize(NPt);
113 ts.resize(NPt);
114 for(int i = 0; i < NPt; ++i) {
115 pts[i] = curveGetPoint(c, i);
116 ts[i] = i / (double)(NPt - 1);
117 }
118 return;
119 }
120 case MSH_SEGM_BEZIER: {
121 int back = -1;
122 int NbCurves = (List_Nbr(c->Control_Points) - 1) / 3;
123 for(int iCurve = 0; iCurve < NbCurves; ++iCurve) {
124 double t1 = (iCurve) / (double)(NbCurves);
125 double t2 = (iCurve + 1) / (double)(NbCurves);
126 SPoint3 pt[4];
127 for(int i = 0; i < 4; i++) { pt[i] = curveGetPoint(c, iCurve * 3 + i); }
128 back = discretizeBezier(tol, discrete, back, pt, t1, t2, iCurve == 0);
129 }
130 break;
131 }
132 case MSH_SEGM_BSPLN: {
133 int back = -1;
134 bool periodic = (c->end == c->beg);
135 int NbControlPoints = List_Nbr(c->Control_Points);
136 int NbCurves = NbControlPoints + (periodic ? -1 : 1);
137 SPoint3 pt[4];
138 for(int iCurve = 0; iCurve < NbCurves; ++iCurve) {
139 double t1 = (iCurve) / (double)(NbCurves);
140 double t2 = (iCurve + 1) / (double)(NbCurves);
141 for(int i = 0; i < 4; i++) {
142 int k;
143 if(periodic) {
144 k = (iCurve - 1 + i) % (NbControlPoints - 1);
145 if(k < 0) k += NbControlPoints - 1;
146 }
147 else {
148 k = std::max(0, std::min(iCurve - 2 + i, NbControlPoints - 1));
149 }
150 pt[i] = curveGetPoint(c, k);
151 }
152 back = discretizeBSpline(tol, discrete, back, pt, t1, t2, iCurve == 0);
153 }
154 break;
155 }
156 case MSH_SEGM_SPLN: {
157 int NbCurves = List_Nbr(c->Control_Points) - 1;
158 SPoint3 pt[4];
159 int back = -1;
160 for(int iCurve = 0; iCurve < NbCurves; ++iCurve) {
161 double t1 = (iCurve) / (double)(NbCurves);
162 double t2 = (iCurve + 1) / (double)(NbCurves);
163 pt[1] = curveGetPoint(c, iCurve);
164 pt[2] = curveGetPoint(c, iCurve + 1);
165 if(iCurve == 0) {
166 if(c->beg == c->end)
167 pt[0] = curveGetPoint(c, NbCurves - 1);
168 else
169 pt[0] = SPoint3(pt[1] * 2 - pt[2]);
170 }
171 else
172 pt[0] = curveGetPoint(c, iCurve - 1);
173 if(iCurve == NbCurves - 1) {
174 if(c->beg == c->end)
175 pt[3] = curveGetPoint(c, 1);
176 else
177 pt[3] = SPoint3(2 * pt[2] - pt[1]);
178 }
179 else
180 pt[3] = curveGetPoint(c, iCurve + 2);
181 back = discretizeCatmullRom(tol, discrete, back, pt, t1, t2, iCurve == 0);
182 }
183 break;
184 }
185 default:
186 Msg::Error("discretizeCurve not implemented for curve type %d", c->Typ);
187 }
188 discrete.sort(pts, ts);
189 }
190
discretize(double tol,std::vector<SPoint3> & dpts,std::vector<double> & ts)191 void gmshEdge::discretize(double tol, std::vector<SPoint3> &dpts,
192 std::vector<double> &ts)
193 {
194 discretizeCurve(c, tol, dpts, ts);
195 }
196