solvespace · Sep 11, 2019 · Sep 11, 2019
Showing with 113 additions and 58 deletions.

+4 −0 CHANGELOG.md

+64 −0 src/mesh.cpp

+6 −0 src/polygon.cpp

+3 −0 src/polygon.h

+36 −58 src/solvespace.cpp
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -70,6 +70,10 @@ 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 calculating area, and faces are selected, calculate area of those faces
+    instead of the closed contour in the sketch.
   * When selecting a point and a line, projected distance to current
     workplane is displayed.
 

diff --git a/src/mesh.cpp b/src/mesh.cpp
@@ -1128,3 +1128,67 @@ 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;
+}
+
+double SMesh::CalculateSurfaceArea(const std::vector<uint32_t> &faces) const {
+    double area = 0.0;
+    for(uint32_t f : faces) {
+        for(const STriangle &t : l) {
+            if(f != t.meta.face) continue;
+            area += t.Area();
+        }
+    }
+    return area;
+}
diff --git a/src/polygon.cpp b/src/polygon.cpp
@@ -87,6 +87,12 @@ double STriangle::SignedVolume() const {
     return a.Dot(b.Cross(c)) / 6.0;
 }
 
+double STriangle::Area() const {
+    Vector ab = a.Minus(b);
+    Vector cb = c.Minus(b);
+    return ab.Cross(cb).Magnitude() / 2.0;
+}
+
 bool STriangle::IsDegenerate() const {
     return a.OnLineSegment(b, c) ||
            b.OnLineSegment(a, c) ||

diff --git a/src/polygon.h b/src/polygon.h
@@ -191,6 +191,7 @@ class STriangle {
     bool Raytrace(const Vector &rayPoint, const Vector &rayDir,
                   double *t, Vector *inters) const;
     double SignedVolume() const;
+    double Area() const;
     bool IsDegenerate() const;
 };
 
@@ -280,6 +281,8 @@ class SMesh {
 
     void PrecomputeTransparency();
     void RemoveDegenerateTriangles();
+    double CalculateVolume() const;
+    double CalculateSurfaceArea(const std::vector<uint32_t> &faces) const;
 
     bool IsEmpty() const;
     void RemapFaces(Group *g, int remap);

diff --git a/src/solvespace.cpp b/src/solvespace.cpp
@@ -772,70 +772,48 @@ void SolveSpaceUI::MenuAnalyze(Command id) {
         }
 
         case Command::VOLUME: {
-            SMesh *m = &(SK.GetGroup(SS.GW.activeGroup)->displayMesh);
-
-            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;
+            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());
+
+            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;
         }
 
         case Command::AREA: {
             Group *g = SK.GetGroup(SS.GW.activeGroup);
+            SS.GW.GroupSelection();
+            auto const &gs = SS.GW.gs;
+            double scale = SS.MmPerUnit();
+
+            if(gs.faces > 0) {
+                std::vector<uint32_t> faces;
+                faces.push_back(gs.face[0].v);
+                if(gs.faces > 1) faces.push_back(gs.face[1].v);
+                double area = g->displayMesh.CalculateSurfaceArea(faces);
+                Message(_("The surface area of the selected faces is:\n\n"
+                          "    %s\n\n"
+                          "Curves have been approximated as piecewise linear.\n"
+                          "This introduces error, typically of around 1%%."),
+                    SS.MmToStringSI(area, /*dim=*/2).c_str());
+                break;
+            }
+
             if(g->polyError.how != PolyError::GOOD) {
                 Error(_("This group does not contain a correctly-formed "
                         "2d closed area. It is open, not coplanar, or self-"