From 915f55aabc201b79abe5f3b70ab9b95ce52a6749 Mon Sep 17 00:00:00 2001 From: EvilSpirit Date: Wed, 11 Sep 2019 10:28:15 +0000 Subject: [PATCH] Show volume of current group alongside total volume. --- CHANGELOG.md | 2 ++ src/mesh.cpp | 53 ++++++++++++++++++++++++++++++++ src/polygon.h | 1 + src/solvespace.cpp | 75 +++++++++++----------------------------------- 4 files changed, 74 insertions(+), 57 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index de34e97e..5dbdf355 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -70,6 +70,8 @@ New measurement/analysis features: * New command for measuring center of mass, with live updates as the sketch changes, "Analyze → Center of Mass". * New option for displaying areas of closed contours. + * When calculating volume of the mesh, volume of the solid from the current + group is now shown alongside total volume of all solids. * When selecting a point and a line, projected distance to current workplane is displayed. diff --git a/src/mesh.cpp b/src/mesh.cpp index 13e0809e..28484204 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1128,3 +1128,56 @@ void SMesh::RemoveDegenerateTriangles() { } l.RemoveTagged(); } + +double SMesh::CalculateVolume() const { + double vol = 0; + for(STriangle tr : l) { + // Translate to place vertex A at (x, y, 0) + Vector trans = Vector::From(tr.a.x, tr.a.y, 0); + tr.a = (tr.a).Minus(trans); + tr.b = (tr.b).Minus(trans); + tr.c = (tr.c).Minus(trans); + + // Rotate to place vertex B on the y-axis. Depending on + // whether the triangle is CW or CCW, C is either to the + // right or to the left of the y-axis. This handles the + // sign of our normal. + Vector u = Vector::From(-tr.b.y, tr.b.x, 0); + u = u.WithMagnitude(1); + Vector v = Vector::From(tr.b.x, tr.b.y, 0); + v = v.WithMagnitude(1); + Vector n = Vector::From(0, 0, 1); + + tr.a = (tr.a).DotInToCsys(u, v, n); + tr.b = (tr.b).DotInToCsys(u, v, n); + tr.c = (tr.c).DotInToCsys(u, v, n); + + n = tr.Normal().WithMagnitude(1); + + // Triangles on edge don't contribute + if(fabs(n.z) < LENGTH_EPS) continue; + + // The plane has equation p dot n = a dot n + double d = (tr.a).Dot(n); + // nx*x + ny*y + nz*z = d + // nz*z = d - nx*x - ny*y + double A = -n.x/n.z, B = -n.y/n.z, C = d/n.z; + + double mac = tr.c.y/tr.c.x, mbc = (tr.c.y - tr.b.y)/tr.c.x; + double xc = tr.c.x, yb = tr.b.y; + + // I asked Maple for + // int(int(A*x + B*y +C, y=mac*x..(mbc*x + yb)), x=0..xc); + double integral = + (1.0/3)*( + A*(mbc-mac)+ + (1.0/2)*B*(mbc*mbc-mac*mac) + )*(xc*xc*xc)+ + (1.0/2)*(A*yb+B*yb*mbc+C*(mbc-mac))*xc*xc+ + C*yb*xc+ + (1.0/2)*B*yb*yb*xc; + + vol += integral; + } + return vol; +} diff --git a/src/polygon.h b/src/polygon.h index c340e136..eed7b155 100644 --- a/src/polygon.h +++ b/src/polygon.h @@ -280,6 +280,7 @@ public: void PrecomputeTransparency(); void RemoveDegenerateTriangles(); + double CalculateVolume() const; bool IsEmpty() const; void RemapFaces(Group *g, int remap); diff --git a/src/solvespace.cpp b/src/solvespace.cpp index b50d8997..0407a941 100644 --- a/src/solvespace.cpp +++ b/src/solvespace.cpp @@ -772,65 +772,26 @@ void SolveSpaceUI::MenuAnalyze(Command id) { } case Command::VOLUME: { - SMesh *m = &(SK.GetGroup(SS.GW.activeGroup)->displayMesh); + Group *g = SK.GetGroup(SS.GW.activeGroup); + double totalVol = g->displayMesh.CalculateVolume(); + std::string msg = ssprintf( + _("The volume of the solid model is:\n\n" + " %s"), + SS.MmToStringSI(totalVol, /*dim=*/3).c_str()); - double vol = 0; - int i; - for(i = 0; i < m->l.n; i++) { - STriangle tr = m->l[i]; - - // Translate to place vertex A at (x, y, 0) - Vector trans = Vector::From(tr.a.x, tr.a.y, 0); - tr.a = (tr.a).Minus(trans); - tr.b = (tr.b).Minus(trans); - tr.c = (tr.c).Minus(trans); - - // Rotate to place vertex B on the y-axis. Depending on - // whether the triangle is CW or CCW, C is either to the - // right or to the left of the y-axis. This handles the - // sign of our normal. - Vector u = Vector::From(-tr.b.y, tr.b.x, 0); - u = u.WithMagnitude(1); - Vector v = Vector::From(tr.b.x, tr.b.y, 0); - v = v.WithMagnitude(1); - Vector n = Vector::From(0, 0, 1); - - tr.a = (tr.a).DotInToCsys(u, v, n); - tr.b = (tr.b).DotInToCsys(u, v, n); - tr.c = (tr.c).DotInToCsys(u, v, n); - - n = tr.Normal().WithMagnitude(1); - - // Triangles on edge don't contribute - if(fabs(n.z) < LENGTH_EPS) continue; - - // The plane has equation p dot n = a dot n - double d = (tr.a).Dot(n); - // nx*x + ny*y + nz*z = d - // nz*z = d - nx*x - ny*y - double A = -n.x/n.z, B = -n.y/n.z, C = d/n.z; - - double mac = tr.c.y/tr.c.x, mbc = (tr.c.y - tr.b.y)/tr.c.x; - double xc = tr.c.x, yb = tr.b.y; - - // I asked Maple for - // int(int(A*x + B*y +C, y=mac*x..(mbc*x + yb)), x=0..xc); - double integral = - (1.0/3)*( - A*(mbc-mac)+ - (1.0/2)*B*(mbc*mbc-mac*mac) - )*(xc*xc*xc)+ - (1.0/2)*(A*yb+B*yb*mbc+C*(mbc-mac))*xc*xc+ - C*yb*xc+ - (1.0/2)*B*yb*yb*xc; - - vol += integral; + SMesh curMesh = {}; + g->thisShell.TriangulateInto(&curMesh); + double curVol = curMesh.CalculateVolume(); + if(curVol > 0.0) { + msg += ssprintf( + _("\nThe volume of current group mesh is:\n\n" + " %s"), + SS.MmToStringSI(curVol, /*dim=*/3).c_str()); } - Message(_("The volume of the solid model is:\n\n" - " %s\n\n" - "Curved surfaces have been approximated as triangles.\n" - "This introduces error, typically of around 1%%."), - SS.MmToStringSI(vol, /*dim=*/3).c_str()); + + msg += _("\n\nCurved surfaces have been approximated as triangles.\n" + "This introduces error, typically of around 1%."); + Message("%s", msg.c_str()); break; }