2

I'm using matlab to calculate the the average velocity in the cross section of a pipe by integrating discrete velocity points over a surface. The points are scattered in a random pattern that form a circle (almost).

I used scatteredinterpolant to create a function relating x and y to v (velocity) in order to create a grid of interpolated values.

F = scatteredInterpolant(x, y, v,'linear');

vq = F(xq,yq); % where xq and yq are a set of query points

The problem I am now having is trying to calculate the the surface area of this function, but only in this circular portion that contains all the scatter points.

The first way I went about this was using the quad2d function.

int = quad2d(@(x,y) F(x,y), min(x), max(x), min(y), max(y), 'MaxFunEvals', 100000);

However this gives is incorrect as it takes the area over a rectangle and a circle.

Now I can probably define the surface area with a circle but in the future I will have to work with more complex shapes so I want to use points that define the boundary of the scatter points.

I'm doing this through triangulation, using the following command.

DT = delaunayTriangulation(x,y);

However I have no idea how I can incorporate these points into a quad2d function. I was hoping someone might have a suggestion or possibly another method I could use in calculating the area over these complex surfaces.

Thanks!

knedlsepp
  • 6,065
  • 3
  • 20
  • 41
user5006725
  • 23
  • 2

1 Answers1

2

You could assume your function to be piecewise linear on your integration area and then integrate it using a midpoint quadrature rule:

For every triangle you compute the midpoint value as the mean of the nodal values and multiply by the triangle's area. You sum it all up to get your integral.

function int = integrate(T, values)
    meanOnTriangle = mean(values(T.ConnectivityList),2);
    int = sum(getElementAreas(T).*meanOnTriangle);        
end

function areas = getElementAreas(T)
    X = @(d) T.Points(T.ConnectivityList(:,d),:);
    d21 = X(2)-X(1);
    d31 = X(3)-X(1);
    areas = abs(1/2*(d21(:,1).*d31(:,2)-d21(:,2).*d31(:,1)));
end

As your goal is the average velocity, you want to compute the following quantity:

averageVelocity = integrate(DT,v)/sum(getElementAreas(DT));
knedlsepp
  • 6,065
  • 3
  • 20
  • 41
  • I was wondering if I had a much more complicated shape such as a heart. Is it possible for triangulation to get the points around the perimeter where the indentation is occurring instead of skipping over it as it is doing currently Or do I need to use a more complicated method. – user5006725 Jun 16 '15 at 17:17
  • @user5006725: That's the problem of triangulating your domain. A delaunay-triangulation will always get you the convex hull (no indentations). You could try [alpha shapes](http://www.mathworks.com/matlabcentral/fileexchange/28851-alpha-shapes) or manually define the boundary and use [constrained delaunay triangulations](http://de.mathworks.com/help/matlab/math/delaunay-triangulation.html#bspqi8a-12). – knedlsepp Jun 16 '15 at 21:30
  • I've created a MATLAB chat room for us to discuss things MATLAB related, or for discussions that span beyond the limitations of a single comment. Visit us when you have time! - http://chat.stackoverflow.com/rooms/81987/matlab – rayryeng Jun 30 '15 at 18:21