@@ -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.
 
 
@@ -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;
+}
@@ -280,6 +280,7 @@ class SMesh {
 
     void PrecomputeTransparency();
     void RemoveDegenerateTriangles();
+    double CalculateVolume() const;
 
     bool IsEmpty() const;
     void RemapFaces(Group *g, int remap);
 
@@ -772,65 +772,26 @@ 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;
         }