How to integrate discrete values over a surface

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]

Do you have a picture? Please attach the data

[Valeria Marangon]

Hi, I'm working on something like this.

Answer 1

[John D'Errico]

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:

https://numfactory.upc.edu/web/Calculo2/P2_Integracio/html/Integral2DQuad.html

[Valeria Marangon]

Hi John, thank you for your answer. That's exactly what I want to do but in my case I know z(x,y) in the center of every quadrangular element and not on every node. That's why I discarded the idea of using a triangular mesh.

I could find a way to calculate node values from element center values but I was hoping to find an easier solution.

[John D'Errico]

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.

quadareas = repmat(diff(rnodes).*(2*rcenternodes*sin(2*pi/(nth-1)/2)),nth-1,1);

Vest1 = sum(quadareas.*zofxy(xcent,ycent),'all')
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]

I've adapted your code to my case and it works perfectly. Thank you so much for your help, John!!

[John D'Errico]

:)

[Michael Obermeyer]

Hi John, pretty late reply on this thread but I encountered a rather similar issue recently [z(x,y) known at the nodes]. A Gaussian quadrature as suggested by you or a mapping of the values to a uniform grid with a successive use of trapz [trapz(y,trapz(x,***))] seem to be the only two viable options.

Your comments on this seem to be correct at first glance except for

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.

Not being a mathematician but only an engineer, I cannot provide an off-the-shelf proof, but please consider the following:

• Triangle {(0,0);(1,0);(1,1)}
• Linear function z(x,y)=1-1*x

Analytical solution:

Your suggested solution:

Please right me if I'm wrong. And no offense - big shout-out for your outstandingly exhaustive answers!

Answer 2

[darova]

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]

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]

• But that is just a complete conjecture at this point.

good

[John D'Errico]

Actually, I was kind of close, though the area of a quadrilateral does figure into the solution. They were not rectangles though.

Answer 3

[Deepanjan Das]

Hi, I am simulating a very simple 2d electrostatic model in Matlab. I have calculated the electric field over a rectangular surface. Now I want to evaluate. So I know the electric field value at each node but not the function. How can I do the integration? Please let me know. Thanks in advance.

[darova]

Are and always the same?

[Deepanjan Das]

No, it depends on the created mesh. Around the boundaries the elements are much smaller in size.