zcy
/
solvespace
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Make surfaces of revolution with control points on the axis 

triangulate correctly; don't screw up generating them, and make
sure that the ratpoly stuff doesn't blow up near the singularity.

[git-p4: depot-paths = "//depot/solvespace/": change = 1953]
solver
Jonathan Westhues 2009-05-18 00:18:32 -08:00
parent 40ed1b7ac1
commit 1d88f96e13
3 changed files with 105 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

23
srf/ratpoly.cpp
View File

@ -301,17 +301,26 @@ Vector SSurface::NormalAt(double u, double v) {
}
void SSurface::ClosestPointTo(Vector p, double *u, double *v, bool converge) {
int i, j;
// A few special cases first; when control points are coincident the
// derivative goes to zero at the conrol points, and would result in
// nonconvergence. We avoid that here, and also guarantee a consistent
// (u, v) (of the infinitely many possible in one parameter).
if(p.Equals(ctrl[0] [0] )) { *u = 0; *v = 0; return; }
if(p.Equals(ctrl[degm][0] )) { *u = 1; *v = 0; return; }
if(p.Equals(ctrl[degm][degn])) { *u = 1; *v = 1; return; }
if(p.Equals(ctrl[0] [degn])) { *u = 0; *v = 1; return; }
// Search for a reasonable initial guess
int i, j;
if(degm == 1 && degn == 1) {
*u = *v = 0; // a plane, perfect no matter what the initial guess
} else {
double minDist = VERY_POSITIVE;
double res = (max(degm, degn) == 2) ? 7.0 : 20.0;
for(i = 0; i < (int)res; i++) {
for(j = 0; j <= (int)res; j++) {
double tryu = (i/res), tryv = (j/res);
int res = (max(degm, degn) == 2) ? 7 : 20;
for(i = 0; i < res; i++) {
for(j = 0; j < res; j++) {
double tryu = (i + 0.5)/res, tryv = (j + 0.5)/res;
Vector tryp = PointAt(tryu, tryv);
double d = (tryp.Minus(p)).Magnitude();
if(d < minDist) {
@ -323,7 +332,7 @@ void SSurface::ClosestPointTo(Vector p, double *u, double *v, bool converge) {
}
}
// Initial guess is in u, v
// Initial guess is in u, v; refine by Newton iteration.
Vector p0;
for(i = 0; i < (converge ? 15 : 3); i++) {
p0 = PointAt(*u, *v);

133
srf/surface.cpp
View File

@ -80,17 +80,25 @@ SSurface SSurface::FromRevolutionOf(SBezier *sb, Vector pt, Vector axis,
Vector ps = p.RotatedAbout(pt, axis, thetas),
pf = p.RotatedAbout(pt, axis, thetaf);
Vector c = ps.ClosestPointOnLine(pt, axis);
Vector ct;
if(ps.Equals(pf)) {
// Degenerate case: a control point lies on the axis of revolution,
// so we get three coincident control points.
ct = ps;
} else {
// Normal case, the control point sweeps out a circle.
Vector c = ps.ClosestPointOnLine(pt, axis);
Vector rs = ps.Minus(c),
rf = pf.Minus(c);
Vector rs = ps.Minus(c),
rf = pf.Minus(c);
Vector ts = axis.Cross(rs),
tf = axis.Cross(rf);
Vector ts = axis.Cross(rs),
tf = axis.Cross(rf);
Vector ct = Vector::AtIntersectionOfLines(ps, ps.Plus(ts),
pf, pf.Plus(tf),
NULL, NULL, NULL);
ct = Vector::AtIntersectionOfLines(ps, ps.Plus(ts),
pf, pf.Plus(tf),
NULL, NULL, NULL);
}
ret.ctrl[i][0] = ps;
ret.ctrl[i][1] = ct;
@ -489,20 +497,36 @@ void SShell::MakeFromRevolutionOf(SBezierLoopSet *sbls, Vector pt, Vector axis,
int color)
{
ZERO(this);
SBezierLoop *sbl;
// Normalize the axis direction so that the direction of revolution
// ends up parallel to the normal of the sketch, on the side of the
// axis where the sketch is.
Vector pto = sbls->point,
ptc = pto.ClosestPointOnLine(pt, axis),
Vector pto;
double md = VERY_NEGATIVE;
for(sbl = sbls->l.First(); sbl; sbl = sbls->l.NextAfter(sbl)) {
SBezier *sb;
for(sb = sbl->l.First(); sb; sb = sbl->l.NextAfter(sb)) {
// Choose the point farthest from the axis; we'll get garbage
// if we choose a point that lies on the axis, for example.
// (And our surface will be self-intersecting if the sketch
// spans the axis, so don't worry about that.)
Vector p = sb->Start();
double d = p.DistanceToLine(pt, axis);
if(d > md) {
md = d;
pto = p;
}
}
}
Vector ptc = pto.ClosestPointOnLine(pt, axis),
up = (pto.Minus(ptc)).WithMagnitude(1),
vp = (sbls->normal).Cross(up);
if(vp.Dot(axis) < 0) {
axis = axis.ScaledBy(-1);
}
SBezierLoop *sbl;
// Now we actually build and trim the surfaces.
for(sbl = sbls->l.First(); sbl; sbl = sbls->l.NextAfter(sbl)) {
int i, j;
SBezier *sb, *prev;
@ -516,11 +540,20 @@ void SShell::MakeFromRevolutionOf(SBezierLoopSet *sbls, Vector pt, Vector axis,
for(sb = sbl->l.First(); sb; sb = sbl->l.NextAfter(sb)) {
Revolved revs;
for(j = 0; j < 4; j++) {
SSurface ss = SSurface::FromRevolutionOf(sb, pt, axis,
(PI/2)*j,
(PI/2)*(j+1));
ss.color = color;
revs.d[j] = surface.AddAndAssignId(&ss);
if(sb->deg == 1 &&
(sb->ctrl[0]).DistanceToLine(pt, axis) < LENGTH_EPS &&
(sb->ctrl[1]).DistanceToLine(pt, axis) < LENGTH_EPS)
{
// This is a line on the axis of revolution; it does
// not contribute a surface.
revs.d[j].v = 0;
} else {
SSurface ss = SSurface::FromRevolutionOf(sb, pt, axis,
(PI/2)*j,
(PI/2)*(j+1));
ss.color = color;
revs.d[j] = surface.AddAndAssignId(&ss);
}
}
hsl.Add(&revs);
}
@ -533,44 +566,54 @@ void SShell::MakeFromRevolutionOf(SBezierLoopSet *sbls, Vector pt, Vector axis,
prev = &(sbl->l.elem[WRAP(i-1, sbl->l.n)]);
for(j = 0; j < 4; j++) {
SCurve sc;
Quaternion qs = Quaternion::From(axis, (PI/2)*j);
// we want Q*(x - p) + p = Q*x + (p - Q*p)
Vector ts = pt.Minus(qs.Rotate(pt));
SCurve sc;
ZERO(&sc);
sc.isExact = true;
sc.exact = sb->TransformedBy(ts, qs);
(sc.exact).MakePwlInto(&(sc.pts));
sc.surfA = revs.d[j];
sc.surfB = revs.d[WRAP(j-1, 4)];
// If this input curve generate a surface, then trim that
// surface with the rotated version of the input curve.
if(revs.d[j].v) {
ZERO(&sc);
sc.isExact = true;
sc.exact = sb->TransformedBy(ts, qs);
(sc.exact).MakePwlInto(&(sc.pts));
sc.surfA = revs.d[j];
sc.surfB = revs.d[WRAP(j-1, 4)];
hSCurve hcb = curve.AddAndAssignId(&sc);
hSCurve hcb = curve.AddAndAssignId(&sc);
STrimBy stb;
stb = STrimBy::EntireCurve(this, hcb, true);
(surface.FindById(sc.surfA))->trim.Add(&stb);
stb = STrimBy::EntireCurve(this, hcb, false);
(surface.FindById(sc.surfB))->trim.Add(&stb);
STrimBy stb;
stb = STrimBy::EntireCurve(this, hcb, true);
(surface.FindById(sc.surfA))->trim.Add(&stb);
stb = STrimBy::EntireCurve(this, hcb, false);
(surface.FindById(sc.surfB))->trim.Add(&stb);
}
SSurface *ss = surface.FindById(sc.surfA);
// And if this input curve and the one after it both generated
// surfaces, then trim both of those by the appropriate
// circle.
if(revs.d[j].v && revsp.d[j].v) {
SSurface *ss = surface.FindById(revs.d[j]);
ZERO(&sc);
sc.isExact = true;
sc.exact = SBezier::From(ss->ctrl[0][0],
ss->ctrl[0][1],
ss->ctrl[0][2]);
sc.exact.weight[1] = ss->weight[0][1];
(sc.exact).MakePwlInto(&(sc.pts));
sc.surfA = revs.d[j];
sc.surfB = revsp.d[j];
ZERO(&sc);
sc.isExact = true;
sc.exact = SBezier::From(ss->ctrl[0][0],
ss->ctrl[0][1],
ss->ctrl[0][2]);
sc.exact.weight[1] = ss->weight[0][1];
(sc.exact).MakePwlInto(&(sc.pts));
sc.surfA = revs.d[j];
sc.surfB = revsp.d[j];
hSCurve hcc = curve.AddAndAssignId(&sc);
hSCurve hcc = curve.AddAndAssignId(&sc);
stb = STrimBy::EntireCurve(this, hcc, false);
(surface.FindById(sc.surfA))->trim.Add(&stb);
stb = STrimBy::EntireCurve(this, hcc, true);
(surface.FindById(sc.surfB))->trim.Add(&stb);
STrimBy stb;
stb = STrimBy::EntireCurve(this, hcc, false);
(surface.FindById(sc.surfA))->trim.Add(&stb);
stb = STrimBy::EntireCurve(this, hcc, true);
(surface.FindById(sc.surfB))->trim.Add(&stb);
}
}
}

2
wishlist.txt
View File

@ -1,7 +1,7 @@
marching algorithm for surface intersection
tangent intersections
surfaces with coincident control points
put back the meshes
assembly
-----