Need surface area of fsurf between two planes

Please help. I have generated a cone with fsurf, and two planes to model a real-world welding problem.

What I am trying to find is the area of the surface of the cone, but only the area between the two planes. I have spent most of the day looking for ways to solve this, but I was not successful. Please help.

conefun = @(x,y,z0,phi) z0 + sqrt(x.^2 + y.^2)*cotd(phi/2);

direction = [0 1 0];
rotate(fsurf(@(x,y) conefun(x,y,0,90)),direction,45)
axis equal
xlim([-1 1.5])
ylim([-1 1])
zlim([ 0 1.5])
xlabel('x-axis')
ylabel('y-axis')
zlabel('z-axis')

hold on

[x y] = meshgrid(-1.5:0.1:1.5);
z = 0.414*x + 0.406;
surf(x,y,z);

[x y] = meshgrid(-1.5:0.1:1.5);
z = 0.414*x + 0.609;
surf(x,y,z);

Thank you very much

Answer 1

[David Goodmanson]

<modified>

Hi Erin,

Here is a totally different method, based on the equation for the cone. It finds the surface area of a cone cut by a plane, then subtracts the area under plane 1 from the area under plane 2. In each case the quantity Sbase is the area of the unperturbed cone, from the lowest height where the plane intersects the cone down to the apex. Phi is the cone opening half-angle, 45 degrees in your case.

I have not downloaded SurfaceIntersection so I can't compare results with the previous answer.

% x = z*tan(phi)*cos(theta);
% y = z*tan(phi)*sin(theta);
% dS = dz/cos(phi) * z*tan(phi)*dtheta     % surface area element
% z = beta*z*tan(phi)*cos(theta) + z0      % the eqn. of the plane
% z1 = z0/(1+beta*tan(phi))                % min and max z where plane cuts cone
% z2 = z0/(1-beta*tan(phi))
% cos(theta) = (z-z0)./(beta*z*tan(phi))   % limits on theta, to stay under the plane

phi = pi/4;
beta = 0.414;
h = [0.406 0.609];

T = tan(phi)/cos(phi);
S = zeros(1,2);
for k = 1:2
  z0 = h(k);
  z1 = z0/(1+beta*tan(phi));
  z2 = z0/(1-beta*tan(phi));
  fun = @(z) 2*z.*acos((z-z0)./(beta*tan(phi)*z));         % modified 
  Sint =  T*integral(fun,z1,z2);
  Sbase = pi*T*z1^2;
  S(k) = Sint + Sbase;
end

Sfinal = diff(S)

[darova]

IS there parameter for rotation of cone?

[Erin Pratt]

Thank you both very much! I choose David's answer as it was very simple to copy and paste into my code, and it works great.

Thanks again

[darova]

How do you know it's correct?

I think my answer is better

[David Goodmanson]

Hello darova,

that's a good question. I like my answer because it gives sensible results and works in a couple of limiting cases. I have not used SurfaceIntersection as of yet. If the answers turn out to be different we should try and determine why that is. If the answers agree, then which one is 'better' is clearly a matter of opinion, nothing more, nothing less.

I didn't use the rotation of the cone because when running Erin's code it didn't seem to make any difference in the 3d plot. It appeared to me that the intent was to have the plane intersect a cone whose axis is vertical, i.e parallel to the z axis, so that's what I did. Anyway If rotation is intended, then a rotation by 45 degrees as specified is going to put the side of the cone on the xy plane, which seems unusual and does not work out well with the planes.

[darova]

I just want to have 2000 reputations points

[Erin Pratt]

darova I tried your routine, but got the error "Undefined function or variable 'SurfaceIntersection'" and so I am not able to try your function. Is there something that I need to download?

[David Goodmanson]

Hello Erin,

In the code I wrote, the axis of the cone is vertical, along the z axis. Did you intend to have any rotation of the cone about the y axis, or was that just for purposes of viewing? And if you did intend for rotation, would that rotation be about the apex of the cone, so the apex would still be in contact with the x-y plane at z =0?

[Erin Pratt]

David

You have it right, the axis of the cone is supposed to be vertical, aligned with the z axis with point of the cone at z=0. When I tried different angles, it seem to yield the correct value, and at an angle of 90 degrees (planes perpendicular to cone axis) the answer was spot on.

Darova

Your method seems very good as well, but the triangularization of the surface area seems to large, the answers are unacceptably off. Is there a way in your code to use more elements, to refine the answer?

I whole-heartedly appreceite both of yours continued involvement. I have to ask what would be the protocol for changing my question a bit? I now realize that the cone in question is not a rotated cone, but rather an eliptical cone with the expression z^2/h^2 = x^2 / a^2 + y^2 / b^2. I was really trying to not take your help for granted, and to figure out how to modify what you have on my own (but I am a Structural Engineer, not a programmer and thus I am a complete newb to Matlab). David, I cant quite figure out how you got the area of the surface in the skewed portion under the planes? I see how you got the total below the lowest point of intersection, but not above between planes. If I understood, then maybe I could modify for my eliptical cone. Darova, I think in your scheme I can just bring in my elipitcal cone, but I need to know how to refine the meshing to get an accurate answer?

Since I owe you reputation points, should I repost a new question with my eliptical cone? That way I can hand out more "Acceptable Answers"?

Please let me know if either of you can further help and how do I thank you within this system?

Thanks!

[David Goodmanson]

Hello Erin,

I found a missing dot in my posted code and have modified the answer accordingly. I don't know how that happened, but it did. The line should be

fun = @(z) 2*z.*acos((z-z0)./(beta*tan(phi)*z));

instead of

fun = @(z) 2*z.*acos((z-z0)/(beta*tan(phi)*z));

(I also made the minor change of renaming a constant).

I wanted to use an independent method for the area (fill the surface of the cone with ten million randomly selected points, find the number of points in between the planes, divide by point surface density). After resinserting the dot, the results agree well. Of course the random point method will always have a bit of variation, at about the third decimal place in this case.

As far as the planes, say you take just the lower of the two planes. That plane cuts the cone between heights z1 (lower) and z2 (upper). In the code, Sint calculates the surface area between the plane and a horizontal surface cutting the cone at height z1. Then Sbase calculates the surface area from height z1 down to the apex, and Sint + Sbase is the area of the cone from the lower plane down to the apex.

Then do the entire process all over again with the upper plane, which has a new z1 and z2, and find the surface area from the upper plane to the apex. Then subtract the two results.

Answer 2

[darova]

Here is what i did

• build all surfaces
• used SurfaceIntersection to get intersection curves (i converted surfaces into patches)

• I got intersection points but in inappropriate order. So i sorted points by angle

• To get the same number of points in each curve i used THIS method. Created patch

Having triangles/faces THIS method was used to calculate surface area

See attached script

[Erin Pratt]

Darova

Thanks again. Yes, I will try your method. Does look most interesting, I will check it out as well.

Thank you so much for putting the time and effort into your response to help me, it was greatly appreceited.

Sincerely

[darova]

I want more reputation points

[Erin Pratt]

I see now. OK, let me check it out.

Sincerely

[David Goodmanson]

Hi Erin, see comment above, in the section for comments to my answer.

[Erin Pratt]

OK, I have the function "SurfaceIntersection()". I am exploring now, thanks!

[darova]

More points

• Since I owe you reputation points, should I repost a new question with my eliptical cone? That way I can hand out more "Acceptable Answers"?

[Erin Pratt]

Darova, yes that helped a bunch. But... Look at the image, there is a missing section. There answer produced was 0.7705 when the correct answer is 0.78098. Any idea on why a slot of area is missing?

[darova]

Sorry, it's because of start/end

Also origin should be specified

rotate(hs,direction,15,[0 0 0])                 % rotate cone

Corrected script. Look now