How to find the equation of the curve that intersects a 3D surface and a plane

Hello, everyone:

I want to calculate the volume of the solid boundary between the surface and two cubes using integral.

I have read this topic 'How do I find the equation of the line that intersects a 3D surface and a plane passing through 3 points of interest on the surface?' and my idea is therefore to fit a polynomial function of the surface (f(x,y)= z = p00 + p10*x + p01*y + p20*x^2 + p11*x*y + p02*y^2)and the intersection between the top plane of the cube and the surface is curveFun: C = p00 + p10*x + p01*y + p20*x^2 + p11*x*y + p02*y^2 and C is the height of the top plane, which is a constant.

This is part of my code:

% to generate the surF, suppose coef are known

zFun = '';
for i=1:np
  for j=1:np
      P = char(vpa(coef(i,j),40));
      if coef(i,j) < 0
          Sign = ' ';
      else
          Sign = ' +';
      end
      I = num2str(i-1);
      J = num2str(j-1);
      zFun = [zFun Sign P '*x^' I '*y^' J];
      %         q = q + 1;
  end
end
sym x,y
surF = matlabFunction(sym(zFun));
% to calculate the volume under the surface
vol = quad2d(surF,xleft,xright,downY,upY);
% to calculate the volume between the surface and top plane 
fun=@(x,y,z)x.^0+y.^0+z.^0;
q = integral3(fun,xleft,xright,Ydown,CurveFun,C,surF);

To calculate the volume under the surface, I can use double integral, which only need the xmin, xmax, ymin, ymax and surfF.

Now for the volume between the surface and the upper cube, I can use the triple integral or I double integral if I know the curve function. Now, how can I construct the curveFun for this triple integral ? The topic mention I can use the Symbolic Math Toolbox to simplify the curve function but I don't know how.

Also, I 'd like to know if there is a better solution to calculate the volume.

I have tried triangulation and the volume under the surface is the sum of these prisms, but the calculation accuracy is largely depended on the number of points I can have.

Many thanks.

Answer 1

[David Sanchez]

Try with the built-in function intersect:

doc intersect

Answer 2

[Josh]

Thanks for you fast reply.

intersect return the common in A and B. However, the surface matrix and the plane may don't have common points.

I draw the intersection curve in the image by:

xRange=(xleft:1:xright);
yRange=(yleft:1:yright);
[X,Y] = meshgrid(xRange,yRange);
zhat = surfF(X,Y);
zdiff = zhat-C;
ContourLine = contours(X, Y, zdiff, [0 0]);

xL = ContourLine (1, 2:end); yL = ContourLine (2, 2:end); zL = interp2(X, Y, zhat, xL, yL); line(xL, yL, zL, 'Color', 'r', 'LineWidth', 3);

I can fit (xL,yL) but I don't want to do that because the fitted function may not match with the "C = p00 + p10*x + p01*y + p20*x^2 + p11*x*y + p02*y^2 "