How to get the cross-sectional area of a PDE model

Now I have a 3D model of a rebar as a .stl file and I want to calculate the cross-sectional area of this rebar at different places. My point is to calculate the area so the contour or slice function doesn't work for me.

I have imported this model to MATLAB as a PDEModel. But when I searched for the cross-sectional area of a PDEModel, I just got how to plot the PDE solution.

Then I used the generateMesh function to get a FEMesh. I think maybe the Surface Intersection function can do it but I don't know how to build a surface from the FEMesh.

So is there any way I can do this?

[darova]

Can you attach the data? Can you show position you want to get crossection?

[Hanshun Yu]

Thanks for your concern. I attached the .stl file here. It's just a cylinder.

About the positoin, I'd like to get cross-sectional area of any section perpendicular to x-axis so it's not so important I think.

Answer 1

[darova]

See this solution

f1 = stlread('Cylinder.stl');

[y,z] = meshgrid(-2:2);             % meshgrid for plane
h = surf(y*0+2,y,z);                % crossection plane
f2 = surf2patch(h,'triangles');     % convert to patch
f2.facecolor= 'r';

patch(f1,'facecolor','y')           % display cylinder
patch(f2)                           % display crossection plane
[~,ff] = SurfaceIntersection(f1,f2);% find crossection line
ff = rmfield(ff,'edges');           % remove edges
patch(ff,'linew',2)                 % display crossection line

alpha(0.3)                          % make surfaces transparent
axis vis3d equal

[Hanshun Yu]

Thank you very much for your answer. I use theSTL File Reader to import the model and get the same result with you. But how should I get the area of cross section? I noticed the 'ff' is a struct, can I get the area directly from it?

[darova]

• extract coordinates from struct

f1 = ff.faces;
yy = ff.vertices(:,2);
zz = ff.vertices(:,3);

You can something like that (data is no ordered)

• order data

[t,r] = cart2pol(yy,zz);        % convert data to polar coordinates
[~,ix] = sort(t);               % sort by angle
plot(yy(ix),zz(ix))

• use polyarea

A = polyarea(yy(ix),zz(ix))
A =
    3.1214

Any idea why the result is not pi=3.1416 ?

[Hanshun Yu]

Thank you very much. I think it can work for now. But what if the surface of cynlinder is not smooth? When I get the model of corroded bars, can I update here to discuss further?

About the answer, I will try the code later so maybe I can get something by then.

[darova]

• But what if the surface of cynlinder is not smooth?

What do you mean?

[Hanshun Yu]

I attached the photo of a corroded rebar. For the model of this rebar, the cross-section is not a circle. Does this effect results of this code?

[darova]

Everything should work ok

[Hanshun Yu]

I see. Thank you again.

[Tim Felix Liesching]

Hey when, I try to run the file Matlab throws an error. My code is a little bit different from yours concerning stlread. stlread() only got the properties Points and ConnectivityList but SurfaceIntersection needs faces and vertices. Do you know how to solve this?

Error:

Unrecognized property 'faces' for class 'triangulation'.

Error in SurfaceIntersection (line 37)

[surface1.faces, surface1.vertices] = freeBoundary(surface1);

Error in untitled_3 (line 55)

[~,ff] = SurfaceIntersection(f1,f2)

>>

f1 = stlread('/Users/tl/Desktop/MATLAB/GUI/STL Dateien/femur.stl');
V = f1.Points;
cl = f1.ConnectivityList;
h1 = patch('Faces',cl,'Vertices',V);

[y,z] =meshgrid(-100:100);
h2 = surf(y,z,z*0+100);
f2 =surf2patch(h2,'triangles')
f2.facecolor = 'r';

patch(f2)
[~,ff] = SurfaceIntersection(f1,f2)
ff = rmfield(ff,'edges')
patch(ff,'linew',2);

[Tim Felix Liesching]

Is it even possible to find the cross-sections, as the cone (see photo) only got vertices at the circular bottom and top plane? There is no vertice data at the shell surface of my stl geometrie. I want to plot several cross sections in 2d depending on where I put my intersection plane.

Best

Tim

[darova]

• My code is a little bit different from yours concerning stlread. stlread() only got the properties Points and ConnectivityList but SurfaceIntersection needs faces and vertices.

try this

f1 = stlread('/Users/tl/Desktop/MATLAB/GUI/STL Dateien/femur.stl');
f1.vertices = f1.points;
f1.faces = f1.connectivitylist;
f1 = rmfield(f1,'points','connectivitylist');

[Tim Felix Liesching]

Thank you very much this works for me!

Do you know why the function is throwing an error when I place the intersection plane at 0 or 50 (lower and upper limit of the cone) ?

h2 = surf(y,z,z*0+50);
h2 = surf(y,z,z*0);

Error:

Error using SurfaceIntersection>IntersectionPolygon (line 709)

Error: in IntersectionPolygon: key is out of bound

Error in SurfaceIntersection>TriangleIntersection2D (line 521)

order = IntersectionPolygon(edgeMat>0, points, dIdx, debug);

Error in SurfaceIntersection (line 263)

[intMatrix(tMsk), intSurface2] = TriangleIntersection2D(V, U, ...

Error in untitled_3 (line 37)

[~,ff] = SurfaceIntersection(str,f2);

This is what my crosse sections look like. Is there a command that fills the cross section with a certain color?

I really appreciate your help.

Best

Tim

[darova]

• Do you know why the function is throwing an error when I place the intersection plane at 0 or 50 (lower and upper limit of the cone) ?

I believe the answer is:

Error: in IntersectionPolygon: key is out of bound

Try

h2 = surf(y,z,z*0+50-0.01);
h2 = surf(y,z,z*0+0.01);

• This is what my crosse sections look like. Is there a command that fills the cross section with a certain color?

Sure, use patch or fill

[Matt]

Hey I'm having a similar issue with this method, I can't find why it doesn't work on my side when trying to "convert" the points/Clist to vertices and faces.

I run the following code :

f1 = stlread('OBJECT.stl');

[y,z] = meshgrid(-20:20);

h = surf(y*0+20,y,z);

f2 = surf2patch(h,'triangles');

f2.facecolor= 'r';

f1.vertices = f1.Points;

f1.faces = f1.ConnectivityList;

f1 = rmfield(f1,'Points','ConnectivityList');

patch(f1,'FaceColor','y')

patch(f2)

but I get this error : Unrecognized property 'vertices' for class 'triangulation'.

and also it seems to Matlab that I don't have enough arguments for patch(f1,'FaceColor','y').

Would appreciate any help !

~Matt

Answer 2

[DGM]

See also this answer using the same FEX tool. This should solve a few of the issues with SurfaceIntersection(). The setup of patch objects is handled internally by the given function, so that solves some usage issues or at least serves as an example. The conversion of unsorted edge lists into sane closed polygons is also handled.

https://www.mathworks.com/matlabcentral/answers/2175386-surface-area-at-given-hight-of-an-stl-file#answer_1562801