How to calculate the volume and centroid of a 3-D geometry that is split by a plane in MATLAB R2023b?

I have a 3-D geometry represented in a STL file, and I need to calculate the volume and centroid of the portion of this geometry which lies below a specified plane. The plane can be at any position and orientation relative to the geometry. How can this be achieved using the Partial Differential Equation Toolbox in MATLAB R2023b?

Answer 1

[MathWorks Support Team]

To calculate the volume and centroid of a portion of a 3-D geometry below a plane in MATLAB R2023b, you can use the Partial Differential Equation Toolbox with the following steps:

1. Import the geometry from the attached STL file using the 'importGeometry' function.
2. Generate a mesh for the imported geometry with the 'generateMesh' function.
3. Compute the incenters of all mesh elements using the 'incenter' method on the mesh object, which will provide the central point of each tetrahedral element.
4. Determine which mesh elements are below the split plane by comparing the z-coordinates of their incenters to the plane equation.
5. Use the 'volume' function in the Partial Differential Equation Toolbox to calculate the total volume of the mesh elements that are identified as being below the plane.
6. Calculate the centroid by averaging the x, y, and z coordinates of the incenters of the mesh elements below the plane.

For higher precision in your calculations, you can reduce the 'HMax' parameter in the 'generateMesh' function. Below is example code for finding the volume and centroid of a sphere split by a plane. Before running this script, download the STL file attached.

%import stl file
model = createpde(1);
importGeometry(model,"sphere.stl");
generateMesh(model,"HMax",0.4);
pdeplot3D(model);

% Obtain connectivity list for triangulation object
mesh = model.Mesh;
pts = mesh.Nodes';
cl = mesh.Elements(1:4,:)';

% Remove unreferenced vertices (only needed for quadratic GeometricOrder)
ind = unique(cl(:));
pp(ind) = 1:length(ind);
cl = pp(cl);
pts = pts(ind,:);

% find centers of triangulation object
tri = triangulation(cl,pts);
surf = triangulation(cl,pts);
centers = incenter(surf)';

% obtain only centers below the plane z = 13
indx = find(centers(3,:)<13);
centers = centers(:,indx);

% find volume and centroid
vol = volume(mesh,indx);
centroid = mean(centers,2);

fprintf("Volume: %f\n",vol);
fprintf("Centroid: [%f, %f, %f]\n",centroid(1),centroid(2),centroid(3));

For more information about these functions from the Partial Differential Equation Toolbox, see their respective documentation pages:

https://www.mathworks.com/help/releases/R2023b/pde/ug/pde.pdemodel.importgeometry.html

https://www.mathworks.com/help/releases/R2023b/pde/ug/pde.pdemodel.generatemesh.html

https://www.mathworks.com/help/releases/R2023b/pde/ug/pde.femesh.volume.html