# How to integrate discrete values over a surface

24 views (last 30 days)
Valeria Marangon on 28 Dec 2019
Commented: John D'Errico on 31 Dec 2019
I'm working on the outputs from a numerical model in which I've used a quadrangular mesh around a cylinder. Each element is defined by its coordinates (x, y) and the corresponding value of the parameter that I want to integrate so I have 3 different arrays containing these information.
My problem is related to the domain because I'm working with two circles, one inside the other (e.g. a plate with a hole) and the data related to the smaller one are NaN. I've tried to use trapz in order to get the difference between the two integrals but of course it didn't work. Is there any other way on how to do this?

darova on 28 Dec 2019
Do you have a picture? Please attach the data
Valeria Marangon on 28 Dec 2019
Hi, I'm working on something like this. John D'Errico on 28 Dec 2019
Edited: John D'Errico on 28 Dec 2019
You are a bit ambiguous here, but I will guess I understand what you are asking.
You have a quad mesh over that domain. At each vertex of the mesh, you have a KNOWN value z(x,y).
Now you want to compute the two dimensional integral of z(x,y) over that domain? So the result will be a single scalar value. If this is your goal, then the answer is pretty simple.
1. Reduce the quadrilateral mesh to a triangle mesh. That just requires you to arbitraily split each quad into a pair of triangles.
2. The integral over that domain is easy, because the integral is just the sum of the integrals over each triangle.
3. The integral over any triangle is easy, because it is just the sum of the z values at each vertex, multiplied by the area of the corresponding triangle, then dividing by 2.
In fact, you can do all of the above in not much more than one line, in a fully vectorized form. I'd probably take a few lines writing it to make it more readable of course.
What I don't know is if this is really your question. (It is often something completely different from what I might guess.) And I don't know for sure how your data is organized. So I won't write code unless I know what the problem really is. Otherwise, I might just be wasting my time.
It might help more to post your data (as a .mat file attachment to a comment.)
Of course, you could also write the solution in a slightly more complex way, working with the quad mesh itself. Then you need to integrate the unknown surface over each quad, then sum the results. Again, not that terribly difficult, IF this is really your goal.
The idea comes up in finite element methods, here is a start to read:

Show 1 older comment
John D'Errico on 29 Dec 2019
That is exactly why I did not spend the time to craft a complete answer. :) When a question is written even slightly confusingly, it always seems to be not what I would expect. At least I was close!
On the other hand, the answer to the modified question is easy. Sort of. In fact the best you can reasonably do is to use a low order Gaussian quadrature over the quad mesh. (Sounds impressive, huh?) The idea is that if you know the value of a function at one point only in a quad mesh, and the point is at the centroid of each quadrilateral, then the integral of z(x,y) over that meshed region is again simply the sum of your estimate over each quad. Basic there. But then, how will you compute the integral of z(x,y) inside a quad element when you know only the center point?
The trick is it is just the value of z(x,y) at the centroid times the area of the quad. In fact, this should be a one point Gaussian quadrature for a quadilateral, so the optimal thing you can do, given only one piece of information per quadrilateral element. It is a LOW order estimate, but since you have only one point per quad, still the best you can do. (Hmm, well, I recall that is the solution for a rect. Not positive if still optimal for a completely general quadrilateral element though.)
Since I don't have your arrays, and I still don't know how you really have things organized in MATLAB, here is the best I can do at the moment:
nr = 11; % radial nodes
nth = 31; % nodes in theta
r1 = 1;
r2 = 2;
rnodes = linspace(r1,r2,nr);
thetanodes = linspace(0,2*pi,nth);
[rmesh,thetamesh] = meshgrid(rnodes,thetanodes);
[xmesh,ymesh] = pol2cart(thetamesh,rmesh);
rcenternodes = rnodes(1:end-1) + diff(rnodes)/2;
thetacenternodes = thetanodes(1:end-1) + diff(thetanodes)/2;
[rcent,thetacent] = meshgrid(rcenternodes,thetacenternodes);
[xcent,ycent] = pol2cart(thetacent,rcent);
figure
% a neat trick here in the plot...
plot(xmesh,ymesh,'k-',xmesh',ymesh','k-',xcent,ycent,'ro')
axis equal
grid on Now I need to fudge up some function for z(x,y) over that region. The nice thing is, I can make it anything I want to compare the final results. Something mildly nonlinear seems a good idea.
zofxy = @(x,y) x.*sin(x+y);
% don't forget that when we work in polar coordinates,
% the differential element is r*dr*dtheta.
zofthetar = @(r,theta) r.*(r.*cos(theta)).*sin(r.*cos(theta) + r.*sin(theta))
format long g
V0 = integral2(zofthetar,r1,r2,0,2*pi)
V0 =
5.36351500290338
So the answer I hope to get is roughly 5.3635... Now let me try it using the quad elements and the low order center point Gaussian quadrature estimate.
First, I need to compute the area of each quad element. Note that they are trapezoidal, with area that varies only as a function of r. The area of a regular trapezoid is easy to compute, just the height of the trap, times the length of the centerline, thus the line of mid-height. Basic geometry there.
So I'll compute the areas of a set of trapezoids that vary along the r dimension, then lazily use repmat to extend that for theta.
Vest1 =
5.35675939530473
Not too bad as an approximation. I'm a bit surprised it is that good, in fact, only 0.12% in error.
(V0 - Vest1)/V0
ans =
0.00125954856003727
Again, the center-point rule is a simple quadrature. If my memory serves me, it will be exact for all functions that are linear over the quad elements, so the error we see will be the error of any deviation from linearity over the small domain of the quad elements. The center-node rule should be a zero'th order Gaussian rule. So it should be exact for polynomials of order 2*0+1 = 1. Late at night for me though, so do some research here to verify my sleepy claim. It is obviously exact for a constant over a fully general quad element.
As long as the quads are sufficiently small, the error will also be small. If you want a true error estimate, it is easy enough to do the numerical analysis to derive an error estimate. Don't trust me as to the claim of error here though. The math is easy, but it is too much the middle of the night for you to trust me.
It should be relatively easy to recraft this for your problem, if I am correct in my guess as to how you are doing things.
Valeria Marangon on 31 Dec 2019
I've adapted your code to my case and it works perfectly. Thank you so much for your help, John!!
John D'Errico on 31 Dec 2019
:)

darova on 28 Dec 2019
If your mesh if fine enough you can just multiply sides  (area of rectangle)
For better precision cross product can be used

John D'Errico on 28 Dec 2019
I think you mistake, that the problem is not to compute the area of a rectangle, but the integral of a function over such a domain, where the function value is given at the vertices of each rectangle in a mesh. But that is just a complete conjecture at this point.
darova on 28 Dec 2019
• But that is just a complete conjecture at this point.
good
John D'Errico on 29 Dec 2019
Actually, I was kind of close, though the area of a quadrilateral does figure into the solution. They were not rectangles though.