Integrating Velocity Over a Complex 2-D Surface

Question

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!

Answer 1

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));