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

Add interpolating splines: both periodic splines (that form a

Browse Source 
loop), and open-ended splines, with their tangents specified at
their endpoints.

Also change constraint solver matrix size to 1024, on the theory
that a power of two will generate better array indexing, and
replace fabs() with my own function that for some reason is
faster.

[git-p4: depot-paths = "//depot/solvespace/": change = 2055]
This commit is contained in:
Jonathan Westhues 2009-10-20 20:46:01 -08:00
parent f6bb680978
commit 2ca5334bdf
18 changed files with 388 additions and 55 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
constraint.cpp
View File

@ -591,12 +591,12 @@ void Constraint::MenuConstrain(int id) {
}
Vector l0 = SK.GetEntity(line->point[0])->PointGetNum(),
l1 = SK.GetEntity(line->point[1])->PointGetNum();
Vector a0 = SK.GetEntity(cubic->point[0])->PointGetNum(),
a3 = SK.GetEntity(cubic->point[3])->PointGetNum();
Vector as = cubic->CubicGetStartNum(),
af = cubic->CubicGetFinishNum();
if(l0.Equals(a0) || l1.Equals(a0)) {
if(l0.Equals(as) || l1.Equals(as)) {
c.other = false;
} else if(l0.Equals(a3) || l1.Equals(a3)) {
} else if(l0.Equals(af) || l1.Equals(af)) {
c.other = true;
} else {
Error("The tangent cubic and line segment must share an "

12
constrainteq.cpp
View File

@ -639,12 +639,12 @@ void ConstraintBase::GenerateReal(IdList<Equation,hEquation> *l) {
EntityBase *cubic = SK.GetEntity(entityA);
EntityBase *line = SK.GetEntity(entityB);
ExprVector endpoint =
SK.GetEntity(cubic->point[other ? 3 : 0])->PointGetExprs();
ExprVector ctrlpoint =
SK.GetEntity(cubic->point[other ? 2 : 1])->PointGetExprs();
ExprVector a = endpoint.Minus(ctrlpoint);
ExprVector a;
if(other) {
a = cubic->CubicGetFinishTangentExprs();
} else {
a = cubic->CubicGetStartTangentExprs();
}
ExprVector b = line->VectorGetExprs();

13
draw.cpp
View File

@ -172,6 +172,7 @@ void GraphicsWindow::GroupSelection(void) {
case Entity::WORKPLANE: (gs.workplanes)++; break;
case Entity::LINE_SEGMENT: (gs.lineSegments)++; break;
case Entity::CUBIC: (gs.cubics)++; break;
// but don't count periodic cubics
case Entity::ARC_OF_CIRCLE:
(gs.circlesOrArcs)++;
@ -203,7 +204,17 @@ void GraphicsWindow::HitTestMakeSelection(Point2d mp) {
for(i = 0; i < SK.entity.n; i++) {
Entity *e = &(SK.entity.elem[i]);
// Don't hover whatever's being dragged.
if(e->h.request().v == pending.point.request().v) continue;
if(e->h.request().v == pending.point.request().v) {
// The one exception is when we're creating a new cubic; we
// want to be able to hover the first point, because that's
// how we turn it into a periodic spline.
if(!e->IsPoint()) continue;
if(!e->h.isFromRequest()) continue;
Request *r = SK.GetRequest(e->h.request());
if(r->type != Request::CUBIC) continue;
if(r->extraPoints < 2) continue;
if(e->h.v != r->h.entity(1).v) continue;
}
d = e->GetDistance(mp);
if(d < 10 && d < dmin) {

4
drawconstraint.cpp
View File

@ -699,8 +699,8 @@ void Constraint::DrawOrGetDistance(Vector *labelPos) {
}
Entity *cubic = SK.GetEntity(entityA);
Vector p =
SK.GetEntity(cubic->point[other ? 3 : 0])->PointGetNum();
Vector p = other ? cubic->CubicGetFinishNum() :
cubic->CubicGetStartNum();
Vector dir = SK.GetEntity(entityB)->VectorGetNum();
Vector out = n.Cross(dir);
textAt = p.Plus(out.WithMagnitude(14/SS.GW.scale));

152
drawentity.cpp
View File

@ -215,6 +215,143 @@ bool Entity::PointIsFromReferences(void) {
return h.request().IsFromReferences();
}
//-----------------------------------------------------------------------------
// Compute a cubic, second derivative continuous, interpolating spline. Same
// routine for periodic splines (in a loop) or open splines (with specified
// end tangents).
//-----------------------------------------------------------------------------
void Entity::ComputeInterpolatingSpline(SBezierList *sbl, bool periodic) {
static const int MAX_N = BandedMatrix::MAX_UNKNOWNS;
int ep = extraPoints;
// The number of unknowns to solve for.
int n = periodic ? 3 + ep : ep;
if(n >= MAX_N) oops();
// The number of on-curve points, one more than the number of segments.
int pts = periodic ? 4 + ep : 2 + ep;
int i, j, a;
// The starting and finishing control points that define our end tangents
// (if the spline isn't periodic), and the on-curve points.
Vector ctrl_s, ctrl_f, pt[MAX_N+4];
if(periodic) {
for(i = 0; i < ep + 3; i++) {
pt[i] = SK.GetEntity(point[i])->PointGetNum();
}
pt[i++] = SK.GetEntity(point[0])->PointGetNum();
} else {
ctrl_s = SK.GetEntity(point[1])->PointGetNum();
ctrl_f = SK.GetEntity(point[ep+2])->PointGetNum();
j = 0;
pt[j++] = SK.GetEntity(point[0])->PointGetNum();
for(i = 2; i <= ep + 1; i++) {
pt[j++] = SK.GetEntity(point[i])->PointGetNum();
}
pt[j++] = SK.GetEntity(point[ep+3])->PointGetNum();
}
// The unknowns that we will be solving for, a set for each coordinate.
double Xx[MAX_N], Xy[MAX_N], Xz[MAX_N];
// For a cubic Bezier section f(t) as t goes from 0 to 1,
// f' (0) = 3*(P1 - P0)
// f' (1) = 3*(P3 - P2)
// f''(0) = 6*(P0 - 2*P1 + P2)
// f''(1) = 6*(P3 - 2*P2 + P1)
for(a = 0; a < 3; a++) {
BandedMatrix bm;
ZERO(&bm);
bm.n = n;
for(i = 0; i < n; i++) {
int im, it, ip;
if(periodic) {
im = WRAP(i - 1, n);
it = i;
ip = WRAP(i + 1, n);
} else {
im = i;
it = i + 1;
ip = i + 2;
}
// All of these are expressed in terms of a constant part, and
// of X[i-1], X[i], and X[i+1]; so let these be the four
// components of that vector;
Vector4 A, B, C, D, E;
// The on-curve interpolated point
C = Vector4::From((pt[it]).Element(a), 0, 0, 0);
// control point one back, C - X[i]
B = C.Plus(Vector4::From(0, 0, -1, 0));
// control point one forward, C + X[i]
D = C.Plus(Vector4::From(0, 0, 1, 0));
// control point two back
if(i == 0 && !periodic) {
A = Vector4::From(ctrl_s.Element(a), 0, 0, 0);
} else {
// pt[im] + X[i-1]
A = Vector4::From(pt[im].Element(a), 1, 0, 0);
}
// control point two forward
if(i == (n - 1) && !periodic) {
E = Vector4::From(ctrl_f.Element(a), 0, 0, 0);
} else {
// pt[ip] - X[i+1]
E = Vector4::From((pt[ip]).Element(a), 0, 0, -1);
}
// Write the second derivatives of each segment, dropping constant
Vector4 fprev_pp = (C.Minus(B.ScaledBy(2))).Plus(A),
fnext_pp = (C.Minus(D.ScaledBy(2))).Plus(E),
eq = fprev_pp.Minus(fnext_pp);
bm.B[i] = -eq.w;
if(periodic) {
bm.A[i][WRAP(i-2, n)] = eq.x;
bm.A[i][WRAP(i-1, n)] = eq.y;
bm.A[i][i] = eq.z;
} else {
// The wrapping would work, except when n = 1 and everything
// wraps to zero...
if(i > 0) bm.A[i][i - 1] = eq.x;
bm.A[i][i] = eq.y;
if(i < (n-1)) bm.A[i][i + 1] = eq.z;
}
}
bm.Solve();
double *X = (a == 0) ? Xx :
(a == 1) ? Xy :
Xz;
memcpy(X, bm.X, n*sizeof(double));
}
for(i = 0; i < pts - 1; i++) {
Vector p0, p1, p2, p3;
if(periodic) {
p0 = pt[i];
int iw = WRAP(i - 1, n);
p1 = p0.Plus(Vector::From(Xx[iw], Xy[iw], Xz[iw]));
} else if(i == 0) {
p0 = pt[0];
p1 = ctrl_s;
} else {
p0 = pt[i];
p1 = p0.Plus(Vector::From(Xx[i-1], Xy[i-1], Xz[i-1]));
}
if(periodic) {
p3 = pt[i+1];
int iw = WRAP(i, n);
p2 = p3.Minus(Vector::From(Xx[iw], Xy[iw], Xz[iw]));
} else if(i == (pts - 2)) {
p3 = pt[pts-1];
p2 = ctrl_f;
} else {
p3 = pt[i+1];
p2 = p3.Minus(Vector::From(Xx[i], Xy[i], Xz[i]));
}
SBezier sb = SBezier::From(p0, p1, p2, p3);
sbl->l.Add(&sb);
}
}
void Entity::GenerateBezierCurves(SBezierList *sbl) {
SBezier sb;
@ -228,15 +365,13 @@ void Entity::GenerateBezierCurves(SBezierList *sbl) {
sbl->l.Add(&sb);
break;
}
case CUBIC: {
Vector p0 = SK.GetEntity(point[0])->PointGetNum();
Vector p1 = SK.GetEntity(point[1])->PointGetNum();
Vector p2 = SK.GetEntity(point[2])->PointGetNum();
Vector p3 = SK.GetEntity(point[3])->PointGetNum();
sb = SBezier::From(p0, p1, p2, p3);
sbl->l.Add(&sb);
case CUBIC:
ComputeInterpolatingSpline(sbl, false);
break;
case CUBIC_PERIODIC:
ComputeInterpolatingSpline(sbl, true);
break;
}
case CIRCLE:
case ARC_OF_CIRCLE: {
@ -476,6 +611,7 @@ void Entity::DrawOrGetDistance(void) {
case CIRCLE:
case ARC_OF_CIRCLE:
case CUBIC:
case CUBIC_PERIODIC:
case TTF_TEXT:
// Nothing but the curve(s).
break;

14
dsc.h
View File

@ -367,4 +367,18 @@ public:
}
};
class BandedMatrix {
public:
static const int MAX_UNKNOWNS = 16;
static const int RIGHT_OF_DIAG = 1;
static const int LEFT_OF_DIAG = 2;
double A[MAX_UNKNOWNS][MAX_UNKNOWNS];
double B[MAX_UNKNOWNS];
double X[MAX_UNKNOWNS];
int n;
void Solve(void);
};
#endif

17
entity.cpp
View File

@ -114,6 +114,23 @@ void EntityBase::ArcGetAngles(double *thetaa, double *thetab, double *dtheta) {
while(*dtheta > (2*PI)) *dtheta -= 2*PI;
}
Vector EntityBase::CubicGetStartNum(void) {
return SK.GetEntity(point[0])->PointGetNum();
}
Vector EntityBase::CubicGetFinishNum(void) {
return SK.GetEntity(point[3+extraPoints])->PointGetNum();
}
ExprVector EntityBase::CubicGetStartTangentExprs(void) {
ExprVector pon = SK.GetEntity(point[0])->PointGetExprs(),
poff = SK.GetEntity(point[1])->PointGetExprs();
return (pon.Minus(poff));
}
ExprVector EntityBase::CubicGetFinishTangentExprs(void) {
ExprVector pon = SK.GetEntity(point[3+extraPoints])->PointGetExprs(),
poff = SK.GetEntity(point[2+extraPoints])->PointGetExprs();
return (pon.Minus(poff));
}
bool EntityBase::IsWorkplane(void) {
return (type == WORKPLANE);
}

10
file.cpp
View File

@ -103,6 +103,7 @@ const SolveSpace::SaveTable SolveSpace::SAVED[] = {
{ 'r', "Request.h.v", 'x', &(SS.sv.r.h.v) },
{ 'r', "Request.type", 'd', &(SS.sv.r.type) },
{ 'r', "Request.extraPoints", 'd', &(SS.sv.r.extraPoints) },
{ 'r', "Request.workplane.v", 'x', &(SS.sv.r.workplane.v) },
{ 'r', "Request.group.v", 'x', &(SS.sv.r.group.v) },
{ 'r', "Request.construction", 'b', &(SS.sv.r.construction) },
@ -120,6 +121,15 @@ const SolveSpace::SaveTable SolveSpace::SAVED[] = {
{ 'e', "Entity.point[1].v", 'x', &(SS.sv.e.point[1].v) },
{ 'e', "Entity.point[2].v", 'x', &(SS.sv.e.point[2].v) },
{ 'e', "Entity.point[3].v", 'x', &(SS.sv.e.point[3].v) },
{ 'e', "Entity.point[4].v", 'x', &(SS.sv.e.point[4].v) },
{ 'e', "Entity.point[5].v", 'x', &(SS.sv.e.point[5].v) },
{ 'e', "Entity.point[6].v", 'x', &(SS.sv.e.point[6].v) },
{ 'e', "Entity.point[7].v", 'x', &(SS.sv.e.point[7].v) },
{ 'e', "Entity.point[8].v", 'x', &(SS.sv.e.point[8].v) },
{ 'e', "Entity.point[9].v", 'x', &(SS.sv.e.point[9].v) },
{ 'e', "Entity.point[10].v", 'x', &(SS.sv.e.point[10].v) },
{ 'e', "Entity.point[11].v", 'x', &(SS.sv.e.point[11].v) },
{ 'e', "Entity.extraPoints", 'd', &(SS.sv.e.extraPoints) },
{ 'e', "Entity.normal.v", 'x', &(SS.sv.e.normal.v) },
{ 'e', "Entity.distance.v", 'x', &(SS.sv.e.distance.v) },
{ 'e', "Entity.workplane.v", 'x', &(SS.sv.e.workplane.v) },

19
group.cpp
View File

@ -636,6 +636,7 @@ void Group::CopyEntity(IdList<Entity,hEntity> *el,
Entity en;
memset(&en, 0, sizeof(en));
en.type = ep->type;
en.extraPoints = ep->extraPoints;
en.h = Remap(ep->h, remap);
en.timesApplied = timesApplied;
en.group = h;
@ -652,12 +653,20 @@ void Group::CopyEntity(IdList<Entity,hEntity> *el,
en.point[1] = Remap(ep->point[1], remap);
break;
case Entity::CUBIC:
en.point[0] = Remap(ep->point[0], remap);
en.point[1] = Remap(ep->point[1], remap);
en.point[2] = Remap(ep->point[2], remap);
en.point[3] = Remap(ep->point[3], remap);
case Entity::CUBIC: {
int i;
for(i = 0; i < 4 + ep->extraPoints; i++) {
en.point[i] = Remap(ep->point[i], remap);
}
break;
}
case Entity::CUBIC_PERIODIC: {
int i;
for(i = 0; i < 3 + ep->extraPoints; i++) {
en.point[i] = Remap(ep->point[i], remap);
}
break;
}
case Entity::CIRCLE:
en.point[0] = Remap(ep->point[0], remap);

4
modify.cpp
View File

@ -310,7 +310,7 @@ hEntity GraphicsWindow::SplitCubic(hEntity he, Vector pinter) {
b01.ClosestPointTo(pinter, &t, true);
b01.SplitAt(t, &b0i, &bi1);
// Add the two line segments this one gets split into.
// Add the two cubic segments this one gets split into.
hRequest r0i = AddRequest(Request::CUBIC, false),
ri1 = AddRequest(Request::CUBIC, false);
// Don't get entities till after adding, realloc issues
@ -343,7 +343,7 @@ hEntity GraphicsWindow::SplitEntity(hEntity he, Vector pinter) {
ret = SplitCircle(he, pinter);
} else if(e->type == Entity::LINE_SEGMENT) {
ret = SplitLine(he, pinter);
} else if(e->type == Entity::CUBIC) {
} else if(e->type == Entity::CUBIC && e->extraPoints == 0) {
ret = SplitCubic(he, pinter);
} else {
Error("Couldn't split this entity; lines, circles, or cubics only.");

92
mouse.cpp
View File

@ -205,12 +205,24 @@ void GraphicsWindow::MouseMoved(double x, double y, bool leftDown,
HitTestMakeSelection(mp);
hRequest hr = pending.point.request();
Vector p0 = SK.GetEntity(hr.entity(1))->PointGetNum();
Vector p3 = SK.GetEntity(hr.entity(4))->PointGetNum();
Vector p1 = p0.ScaledBy(2.0/3).Plus(p3.ScaledBy(1.0/3));
SK.GetEntity(hr.entity(2))->PointForceTo(p1);
Vector p2 = p0.ScaledBy(1.0/3).Plus(p3.ScaledBy(2.0/3));
SK.GetEntity(hr.entity(3))->PointForceTo(p2);
if(pending.point.v == hr.entity(4).v) {
// The very first segment; dragging final point drags both
// tangent points.
Vector p0 = SK.GetEntity(hr.entity(1))->PointGetNum(),
p3 = SK.GetEntity(hr.entity(4))->PointGetNum(),
p1 = p0.ScaledBy(2.0/3).Plus(p3.ScaledBy(1.0/3)),
p2 = p0.ScaledBy(1.0/3).Plus(p3.ScaledBy(2.0/3));
SK.GetEntity(hr.entity(1+1))->PointForceTo(p1);
SK.GetEntity(hr.entity(1+2))->PointForceTo(p2);
} else {
// A subsequent segment; dragging point drags only final
// tangent point.
int i = SK.GetEntity(hr.entity(0))->extraPoints;
Vector pn = SK.GetEntity(hr.entity(4+i))->PointGetNum(),
pnm2 = SK.GetEntity(hr.entity(2+i))->PointGetNum(),
pnm1 = (pn.Plus(pnm2)).ScaledBy(0.5);
SK.GetEntity(hr.entity(3+i))->PointForceTo(pnm1);
}
SS.MarkGroupDirtyByEntity(pending.point);
break;
@ -282,6 +294,7 @@ void GraphicsWindow::MouseMoved(double x, double y, bool leftDown,
void GraphicsWindow::ClearPending(void) {
memset(&pending, 0, sizeof(pending));
SS.later.showTW = true;
}
void GraphicsWindow::MouseMiddleOrRightDown(double x, double y) {
@ -319,15 +332,19 @@ void GraphicsWindow::MouseRightUp(double x, double y) {
if(orig.startedMoving) return;
if(context.active) return;
context.active = true;
if(pending.operation == DRAGGING_NEW_LINE_POINT) {
if(pending.operation == DRAGGING_NEW_LINE_POINT ||
pending.operation == DRAGGING_NEW_CUBIC_POINT)
{
// Special case; use a right click to stop drawing lines, since
// a left click would draw another one. This is quicker and more
// intuitive than hitting escape.
// intuitive than hitting escape. Likewise for new cubic segments.
ClearPending();
return;
}
context.active = true;
GroupSelection();
if(hover.IsEmpty() && gs.n == 0 && gs.constraints == 0) {
// No reason to display a context menu.
@ -591,7 +608,7 @@ void GraphicsWindow::MouseLeftDown(double mx, double my) {
pending.operation = DRAGGING_NEW_LINE_POINT;
pending.point = hr.entity(2);
pending.description = "click to place next point of line";
pending.description = "click next point of line, or press Esc";
SK.GetEntity(pending.point)->PointForceTo(v);
break;
@ -675,7 +692,7 @@ void GraphicsWindow::MouseLeftDown(double mx, double my) {
pending.operation = DRAGGING_NEW_CUBIC_POINT;
pending.point = hr.entity(4);
pending.description = "click to place next point of cubic";
pending.description = "click next point of cubic, or press Esc";
break;
case MNU_WORKPLANE:
@ -734,11 +751,60 @@ void GraphicsWindow::MouseLeftDown(double mx, double my) {
break;
case DRAGGING_NEW_ARC_POINT:
case DRAGGING_NEW_CUBIC_POINT:
ConstrainPointByHovered(pending.point);
ClearPending();
break;
case DRAGGING_NEW_CUBIC_POINT: {
hRequest hr = pending.point.request();
Request *r = SK.GetRequest(hr);
if(hover.entity.v == hr.entity(1).v && r->extraPoints >= 2) {
// They want the endpoints coincident, which means a periodic
// spline instead.
r->type = Request::CUBIC_PERIODIC;
// Remove the off-curve control points, which are no longer
// needed here; so move [2,ep+1] down, skipping first pt.
int i;
for(i = 2; i <= r->extraPoints+1; i++) {
SK.GetEntity(hr.entity((i-1)+1))->PointForceTo(
SK.GetEntity(hr.entity(i+1))->PointGetNum());
}
// and move ep+3 down by two, skipping both
SK.GetEntity(hr.entity((r->extraPoints+1)+1))->PointForceTo(
SK.GetEntity(hr.entity((r->extraPoints+3)+1))->PointGetNum());
r->extraPoints -= 2;
// And we're done.
SS.MarkGroupDirty(r->group);
SS.later.generateAll = true;
ClearPending();
break;
}
if(ConstrainPointByHovered(pending.point)) {
ClearPending();
break;
}
Entity e;
if(r->extraPoints >= arraylen(e.point) - 4) {
ClearPending();
break;
}
(SK.GetRequest(hr)->extraPoints)++;
SS.GenerateAll(-1, -1);
int ep = r->extraPoints;
Vector last = SK.GetEntity(hr.entity(3+ep))->PointGetNum();
SK.GetEntity(hr.entity(2+ep))->PointForceTo(last);
SK.GetEntity(hr.entity(3+ep))->PointForceTo(v);
SK.GetEntity(hr.entity(4+ep))->PointForceTo(v);
pending.point = hr.entity(4+ep);
break;
}
case DRAGGING_NEW_LINE_POINT: {
if(ConstrainPointByHovered(pending.point)) {
ClearPending();
@ -755,7 +821,7 @@ void GraphicsWindow::MouseLeftDown(double mx, double my) {
// And drag an endpoint of the new line segment
pending.operation = DRAGGING_NEW_LINE_POINT;
pending.point = hr.entity(2);
pending.description = "click to place next point of next line";
pending.description = "click next point of line, or press Esc";
break;
}

9
request.cpp
View File

@ -50,7 +50,12 @@ void Request::Generate(IdList<Entity,hEntity> *entity,
case Request::CUBIC:
et = Entity::CUBIC;
points = 4;
points = 4 + extraPoints;
break;
case Request::CUBIC_PERIODIC:
et = Entity::CUBIC_PERIODIC;
points = 3 + extraPoints;
break;
case Request::TTF_TEXT:
@ -64,6 +69,7 @@ void Request::Generate(IdList<Entity,hEntity> *entity,
// Generate the entity that's specific to this request.
e.type = et;
e.extraPoints = extraPoints;
e.group = group;
e.style = style;
e.workplane = workplane;
@ -157,6 +163,7 @@ char *Request::DescriptionString(void) {
case DATUM_POINT: s = "datum-point"; break;
case LINE_SEGMENT: s = "line-segment"; break;
case CUBIC: s = "cubic-bezier"; break;
case CUBIC_PERIODIC: s = "periodic-cubic"; break;
case CIRCLE: s = "circle"; break;
case ARC_OF_CIRCLE: s = "arc-of-circle"; break;
case TTF_TEXT: s = "ttf-text"; break;

12
sketch.h
View File

@ -248,11 +248,13 @@ public:
static const int DATUM_POINT = 101;
static const int LINE_SEGMENT = 200;
static const int CUBIC = 300;
static const int CUBIC_PERIODIC = 301;
static const int CIRCLE = 400;
static const int ARC_OF_CIRCLE = 500;
static const int TTF_TEXT = 600;
int type;
int extraPoints;
hEntity workplane; // or Entity::FREE_IN_3D
hGroup group;
@ -302,6 +304,7 @@ public:
static const int WORKPLANE = 10000;
static const int LINE_SEGMENT = 11000;
static const int CUBIC = 12000;
static const int CUBIC_PERIODIC = 12001;
static const int CIRCLE = 13000;
static const int ARC_OF_CIRCLE = 14000;
static const int TTF_TEXT = 15000;
@ -313,7 +316,8 @@ public:
// When it comes time to draw an entity, we look here to get the
// defining variables.
hEntity point[4];
hEntity point[12];
int extraPoints;
hEntity normal;
hEntity distance;
// The only types that have their own params are points, normals,
@ -386,6 +390,11 @@ public:
ExprVector NormalExprsV(void);
ExprVector NormalExprsN(void);
Vector CubicGetStartNum(void);
Vector CubicGetFinishNum(void);
ExprVector CubicGetStartTangentExprs(void);
ExprVector CubicGetFinishTangentExprs(void);
void AddEq(IdList<Equation,hEquation> *l, Expr *expr, int index);
void GenerateEquations(IdList<Equation,hEquation> *l);
};
@ -423,6 +432,7 @@ public:
bool IsVisible(void);
bool PointIsFromReferences(void);
void ComputeInterpolatingSpline(SBezierList *sbl, bool periodic);
void GenerateBezierCurves(SBezierList *sbl);
void GenerateEdges(SEdgeList *el, bool includingConstruction=false);

7
solvespace.h
View File

@ -33,6 +33,9 @@ inline double WRAP_SYMMETRIC(double v, double n) {
return v;
}
// Why is this faster than the library function?
inline double ffabs(double v) { return (v > 0) ? v : (-v); }
#define SWAP(T, a, b) do { T temp = (a); (a) = (b); (b) = temp; } while(0)
#define ZERO(v) memset((v), 0, sizeof(*(v)))
#define CO(v) (v).x, (v).y, (v).z
@ -234,8 +237,8 @@ bool StringEndsIn(char *str, char *ending);
class System {
public:
static const int MAX_UNKNOWNS = 1000;
static const int MAX_DRAGGED = 4;
static const int MAX_UNKNOWNS = 1024;
static const int MAX_DRAGGED = 4;
EntityList entity;
ParamList param;

12
system.cpp
View File

@ -171,15 +171,15 @@ bool System::SolveLinearSystem(double X[], double A[][MAX_UNKNOWNS],
// greater. First, find a pivot (between rows i and N-1).
max = 0;
for(ip = i; ip < n; ip++) {
if(fabs(A[ip][i]) > max) {
if(ffabs(A[ip][i]) > max) {
imax = ip;
max = fabs(A[ip][i]);
max = ffabs(A[ip][i]);
}
}
// Don't give up on a singular matrix unless it's really bad; the
// assumption code is responsible for identifying that condition,
// so we're not responsible for reporting that error.
if(fabs(max) < 1e-20) return false;
if(ffabs(max) < 1e-20) return false;
// Swap row imax with row i
for(jp = 0; jp < n; jp++) {
@ -201,7 +201,7 @@ bool System::SolveLinearSystem(double X[], double A[][MAX_UNKNOWNS],
// We've put the matrix in upper triangular form, so at this point we
// can solve by back-substitution.
for(i = n - 1; i >= 0; i--) {
if(fabs(A[i][i]) < 1e-20) return false;
if(ffabs(A[i][i]) < 1e-20) return false;
temp = B[i];
for(j = n - 1; j > i; j--) {
@ -294,7 +294,7 @@ bool System::NewtonSolve(int tag) {
if(isnan(mat.B.num[i])) {
return false;
}
if(fabs(mat.B.num[i]) > CONVERGE_TOLERANCE) {
if(ffabs(mat.B.num[i]) > CONVERGE_TOLERANCE) {
converged = false;
break;
}
@ -483,7 +483,7 @@ int System::Solve(Group *g, int *dof, List<hConstraint> *bad,
didnt_converge:
SK.constraint.ClearTags();
for(i = 0; i < eq.n; i++) {
if(fabs(mat.B.num[i]) > CONVERGE_TOLERANCE || isnan(mat.B.num[i])) {
if(ffabs(mat.B.num[i]) > CONVERGE_TOLERANCE || isnan(mat.B.num[i])) {
// This constraint is unsatisfied.
if(!mat.eq[i].isFromConstraint()) continue;

17
textwin.cpp
View File

@ -335,11 +335,22 @@ void TextWindow::DescribeSelection(void) {
SS.MmToString((p1.Minus(p0).Magnitude())));
break;
}
case Entity::CUBIC_PERIODIC:
case Entity::CUBIC:
Printf(false, "%FtCUBIC BEZIER CURVE%E");
for(i = 0; i <= 3; i++) {
int pts;
if(e->type == Entity::CUBIC_PERIODIC) {
Printf(false, "%FtPERIODIC C2 CUBIC SPLINE%E");
pts = (3 + e->extraPoints);
} else if(e->extraPoints > 0) {
Printf(false, "%FtINTERPOLATING C2 CUBIC SPLINE%E");
pts = (4 + e->extraPoints);
} else {
Printf(false, "%FtCUBIC BEZIER CURVE%E");
pts = 4;
}
for(i = 0; i < pts; i++) {
p = SK.GetEntity(e->point[i])->PointGetNum();
Printf((i==0), " p%c = " PT_AS_STR, '0'+i, COSTR(p));
Printf((i==0), " p%d = " PT_AS_STR, i, COSTR(p));
}
break;

40
util.cpp
View File

@ -125,6 +125,46 @@ void MakeMatrix(double *mat, double a11, double a12, double a13, double a14,
mat[15] = a44;
}
//-----------------------------------------------------------------------------
// Solve a mostly banded matrix. In a given row, there are LEFT_OF_DIAG
// elements to the left of the diagonal element, and RIGHT_OF_DIAG elements to
// the right (so that the total band width is LEFT_OF_DIAG + RIGHT_OF_DIAG + 1).
// There also may be elements in the last two columns of any row. We solve
// without pivoting.
//-----------------------------------------------------------------------------
void BandedMatrix::Solve(void) {
int i, ip, j, jp;
double temp;
// Reduce the matrix to upper triangular form.
for(i = 0; i < n; i++) {
for(ip = i+1; ip < n && ip <= (i + LEFT_OF_DIAG); ip++) {
temp = A[ip][i]/A[i][i];
for(jp = i; jp < (n - 2) && jp <= (i + RIGHT_OF_DIAG); jp++) {
A[ip][jp] -= temp*(A[i][jp]);
}
A[ip][n-2] -= temp*(A[i][n-2]);
A[ip][n-1] -= temp*(A[i][n-1]);
B[ip] -= temp*B[i];
}
}
// And back-substitute.
for(i = n - 1; i >= 0; i--) {
temp = B[i];
if(i < n-1) temp -= X[n-1]*A[i][n-1];
if(i < n-2) temp -= X[n-2]*A[i][n-2];
for(j = min(n - 3, i + RIGHT_OF_DIAG); j > i; j--) {
temp -= X[j]*A[i][j];
}
X[i] = temp / A[i][i];
}
}
const Quaternion Quaternion::IDENTITY = { 1, 0, 0, 0 };
Quaternion Quaternion::From(double w, double vx, double vy, double vz) {

1
wishlist.txt
View File

@ -1,7 +1,6 @@
multi-drag
copy and paste
interpolating splines
associative entities from solid model, as a special group
-----