Thread Subject: Area integral over a triangle in MATLAB -- is numerical integration possible?

Hi,

I'm wondering if it is possible to perform a double integral in MATLAB like

\int \int f(x,y) dx dy

where the domain of integration is a triangle, whose vertex coordinates are known?

A full question is posted here: http://www.physicsforums.com/showthread.php?t=388101.

Thanks in advance!

"Vivek Saxena" <maverick280857@yahoo.com> wrote in message
news:ho1k5u$hsc$1@fred.mathworks.com...
> Hi,
>
> I'm wondering if it is possible to perform a double integral in MATLAB
> like
>
> \int \int f(x,y) dx dy
>
> where the domain of integration is a triangle, whose vertex coordinates
> are known?
>
> A full question is posted here:
> http://www.physicsforums.com/showthread.php?t=388101.

Look at QUAD2D.

--
Steve Lord
slord@mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ

"Steven Lord" <slord@mathworks.com> wrote in message <ho1kqh$qvl$1@fred.mathworks.com>...
>
> "Vivek Saxena" <maverick280857@yahoo.com> wrote in message
> news:ho1k5u$hsc$1@fred.mathworks.com...
> > Hi,
> >
> > I'm wondering if it is possible to perform a double integral in MATLAB
> > like
> >
> > \int \int f(x,y) dx dy
> >
> > where the domain of integration is a triangle, whose vertex coordinates
> > are known?
> >
> > A full question is posted here:
> > http://www.physicsforums.com/showthread.php?t=388101.
>
> Look at QUAD2D.
>
> --
> Steve Lord
> slord@mathworks.com
> comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
>

Thanks Steve. The approach I thought to define c(x) and d(x) is based on the equations of the sides of the triangle (defining the boundaries of the area element). But there are several possibilities depending on how the triangle is oriented with respect to the axes: an edge could have a zero slope, a nonzero finite slope or an infinite slope, and the different permutations of the ordering of the vertices will probably also play a role in defining c(x) and d(x) differently, for different triangles.

Is there a neater way to employ QUAD2D to find the integral of f(x,y) over a triangular region with arbitrary vertex locations (x1, y1), (x2, y2) and (x3, y3)?

"Vivek Saxena" <maverick280857@yahoo.com> wrote in message <ho1k5u$hsc$1@fred.mathworks.com>...
> I'm wondering if it is possible to perform a double integral in MATLAB like
> \int \int f(x,y) dx dy
> where the domain of integration is a triangle, whose vertex coordinates are known?

  If it happens you don't have quad2d, you can use the older dblquad which integrates over rectangular regions by making a change of variable from x,y to u,v defined by:

 x = (1-u)*x2 + u*((1-v)*x2 + v*x3)
 y = (1-u)*y2 + u*((1-v)*y2 + v*y3)

where (x1,y1), (x2,y2), and (x3,y3) are the triangular region's three vertices. Of course you must change dx dy in the double integral to du dv with the proper Jacobian factor. What was the triangular region in x,y space then becomes a unit square in u,v space - the variables u and v each range from 0 to 1.

  Of course it is more convenient to use quad2d if you have it.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <ho1ocu$ldr$1@fred.mathworks.com>...
> "Vivek Saxena" <maverick280857@yahoo.com> wrote in message <ho1k5u$hsc$1@fred.mathworks.com>...
> > I'm wondering if it is possible to perform a double integral in MATLAB like
> > \int \int f(x,y) dx dy
> > where the domain of integration is a triangle, whose vertex coordinates are known?
>
> If it happens you don't have quad2d, you can use the older dblquad which integrates over rectangular regions by making a change of variable from x,y to u,v defined by:
>
> x = (1-u)*x2 + u*((1-v)*x2 + v*x3)
> y = (1-u)*y2 + u*((1-v)*y2 + v*y3)
>
> where (x1,y1), (x2,y2), and (x3,y3) are the triangular region's three vertices. Of course you must change dx dy in the double integral to du dv with the proper Jacobian factor. What was the triangular region in x,y space then becomes a unit square in u,v space - the variables u and v each range from 0 to 1.
>
> Of course it is more convenient to use quad2d if you have it.
>
> Roger Stafford

Thanks Roger. But didn't you mean

x = (1-u)*x1 + u*[(1-v)*x2 + v*x3]
y = (1-u)*y1 + u*[(1-v)*y2 + v*y3]

because the coordinate map you gave doesn't involve x1 and y1. But even so, this becomes a triangle in (u, v) space because

u = 0, v = 0 => x = x1, y = y1
u = 0, v = 1 => x = x1, y = y1
u = 1, v = 0 => x = x2, y = y2
u = 1, v = 1 => x = x3, y = y3

Only the last 3 coordinates are meaningful essentially, because if u = 0, the way this map is constructed, the value of v becomes meaningless. Perhaps you meant something else?

But I think this will do, because this transforms the triangle into a right angled triangle in (u, v) space with the vertices (0, 0), (0, 1), and (1,1).

Now I just have to write the integrand as f(x(u,v), y(u,v)), multiply it by the Jacobian factor, and integrate from u = 0 to 1 and v = 0 to u (and not v = 0 to 1). Is it possible to do this symbolically? That is, what we now have to do is

\int_{u=0}^{u=1} \int_{v=0}^{v=u} g(u, v) du dv

What limits does one specify to quad2d here?

"Vivek Saxena" <maverick280857@yahoo.com> wrote in message <ho1pvj$cbr$1@fred.mathworks.com>...
> Thanks Roger. But didn't you mean
>
> x = (1-u)*x1 + u*[(1-v)*x2 + v*x3]
> y = (1-u)*y1 + u*[(1-v)*y2 + v*y3]
>
> because the coordinate map you gave doesn't involve x1 and y1. But even so, this becomes a triangle in (u, v) space because
>
> u = 0, v = 0 => x = x1, y = y1
> u = 0, v = 1 => x = x1, y = y1
> u = 1, v = 0 => x = x2, y = y2
> u = 1, v = 1 => x = x3, y = y3
>
> Only the last 3 coordinates are meaningful essentially, because if u = 0, the way this map is constructed, the value of v becomes meaningless. Perhaps you meant something else?

  You are quite right in your first assertion, Vivek. It was a typing error and should be as you say,

x = (1-u)*x1 + u*((1-v)*x2 + v*x3)
y = (1-u)*y1 + u*((1-v)*y2 + v*y3)

  However I disagree with your second statement. The interior of the triangle in x and y becomes the full interior of the unit square in u & v. At the far end the point (x1,y1) becomes stretched out to the line segment connecting (0,0) and (0,1) in u,v space. Just shy of (x1,y1) a very short strip in x,y space is stretched to the full unit length strip in u,v space. The transformation is valid because at the same time the Jacobian shrinks toward zero in compensation as you approach the point (x1,y1). Of course the inverse transformation from u,v to x,y becomes singular at that point but it doesn't matter in finding the double integral. It only means that dblquad will be a little wasteful in the number of points it uses toward the u = 0 end.

  If you doubt the above statement, try this:

x1 = randn; y1 = randn; % Random vertices for triangle in x,y
x2 = randn; y2 = randn;
x3 = randn; y3 = randn;
N = 2000;
u = rand(1,N); % u & v fill the unit square in u,v space
v = rand(1,N);
x = (1-u)*x1 + u.*((1-v)*x2 + v*x3); % x and y will fill the triangle
y = (1-u)*y1 + u.*((1-v)*y2 + v*y3);
plot([x1,x2,x3,x1],[y1,y2,y3,y1],'r-',x,y,'y.')
axis equal

You will see that the red triangle becomes nicely filled up with yellow dots even though they become tightly packed toward the (x1,y1) vertex, (where the Jacobian is getting smaller.)

Roger Stafford

"Vivek Saxena" <maverick280857@yahoo.com> wrote in message <ho1k5u$hsc$1@fred.mathworks.com>...
> Hi,
>
> I'm wondering if it is possible to perform a double integral in MATLAB like
>
> \int \int f(x,y) dx dy
>
> where the domain of integration is a triangle, whose vertex coordinates are known?
>
> A full question is posted here: http://www.physicsforums.com/showthread.php?t=388101.

It looks like a wave equation. Usually in FEM integration perform using Gaussian quadrature formula. On triangle, it involve barycentric coordinates. There are several schemes more or less complex and accurate.

Very accurate integration such as QUAD2 is not necessary in FEM, but take a lot of time.

Bruno

"Vivek Saxena" <maverick280857@yahoo.com> wrote in message <ho1rq8$8rk$1@fred.mathworks.com>...
> But I think this will do, because this transforms the triangle into a right angled triangle in (u, v) space with the vertices (0, 0), (0, 1), and (1,1).
>
> Now I just have to write the integrand as f(x(u,v), y(u,v)), multiply it by the Jacobian factor, and integrate from u = 0 to 1 and v = 0 to u (and not v = 0 to 1). Is it possible to do this symbolically? That is, what we now have to do is
>
> \int_{u=0}^{u=1} \int_{v=0}^{v=u} g(u, v) du dv
>
> What limits does one specify to quad2d here?

  No Vivek, you should integrate over the full square u = 0 to u = 1 and v = 0 to v = 1. dblquad will not let you do otherwise. It is the full u,v square that the x,y triangle is mapped into by the given transformation. It does not matter that it has a singularity at the (x1,y1) vertex.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <ho1u8b$dfi$1@fred.mathworks.com>...
> No Vivek, you should integrate over the full square u = 0 to u = 1 and v = 0 to v = 1. dblquad will not let you do otherwise. It is the full u,v square that the x,y triangle is mapped into by the given transformation. It does not matter that it has a singularity at the (x1,y1) vertex.
>
> Roger Stafford

Hi Roger, I think I see your point. But what if I use quad2d?

So, if I understand you correctly, the triangle is mapped to the square {(u,v): 0<= u <= 1, 0 <= v <= 1} in (u,v) space.

"Vivek Saxena" <maverick280857@yahoo.com> wrote in message <ho24mk$mb1$1@fred.mathworks.com>...
> Hi Roger, I think I see your point. But what if I use quad2d?

  If you use quad2d, then you could use the original x,y coordinates directly without any transformation, though that is a bit awkward because in general you have to break the integration into two separate ranges which depend on how the three points of the triangle are oriented. The triangle's edges may not be aligned with either of the coordinate axes.

  It might be more convenient to use an affine transformation that transforms your triangle into another triangle that is oriented along the coordinate directions so that separate integration ranges would not be required. For example, you could say

 x = x1 + u*(x2-x1) + v*(x3-x1)
 y = y1 + u*(y2-y1) + v*(y3-y1)

This transforms the triangular region u >= 0, v >= 0, u+v <= 1 into your given x,y triangle. Thus your quad2d in the u,v domain could have simple limits: u from 0 to 1 and v from 0 to 1-u, and that is within quad2d's capabilities.

  Of course you must still use the correct Jacobian factor in the integral if you do this transformation, but it is now just a constant quantity, namely the ratio between the two triangles' respective areas.

> So, if I understand you correctly, the triangle is mapped to the square {(u,v): 0<= u <= 1, 0 <= v <= 1} in (u,v) space.

Yes, that is correct (for the earlier transformation using dblquad.)

Roger Stafford

>
> Yes, that is correct (for the earlier transformation using dblquad.)
>
> Roger Stafford

Hi again Roger. I thought about this a bit, and I'm unable to convince myself that the nonlinear transformation you wrote really does map a triangle into a square. I guess you're right, but when I did the integral, the values all came out wrong (I have an analytical solution to the final problem with which I can compare). And since this is the only additional term in my calculation, I have no reason to suspect the rest of the code which has been working so far.

Anyhow, I ran the following code

clear, close all;
clc;
x1 = 0, x2 = 1, x3 = 1, y1 = 0, y2 = 0, y3 = 1;
c = 1;
for u=0:0.01:1
    for v=0:0.01:1
        x(c) = (1-u)*x1 + u*((1-v)*x2+v*x3);
        y(c) = (1-u)*y1 + u*((1-v)*y2+v*y3);
        c = c + 1;
    end
end
scatter(x,y);

That is I begin with a triangle which is right angled, with its vertices in (x, y) space as (0,0), (0,1) and (1,1) for simplicity. Now I use this transformation to construct x and y points for u and v ranging from 0 to 1. This gives me the coordinates of all the points inside the triangle (as required). But here, I have assumed that the corresponding region in (u, v) space is a square (because I allow u and v to run freely). So the forward transformation is a map from a square to a triangle. I am convinced about that. But how do I prove that the inverse transformation maps the triangle to the square {(u,v) : 0<=u<=1, 0<=v<=1}?

Perhaps this is a trivial question, but I can't seem to get my head around it. I am using dblquad, and my functions look like this:

function f = abcintegrand(u,v,a,b,c,k_0,x1,x2,x3,y1,y2,y3)

J = sqrt(-1);
argexp = (1-u).*x1 + u.*((1-v).*x2 + v.*x3);
jacob = (-x1+(1-v).*x2+v.*x3).*(-u.*y2+u.*y3)-(-y1+(1-v).*y2+v.*y3).*(-u.*x2+u.*x3);
p = (a + b*((1-u).*x1 + u.*((1-v).*x2 + v.*x3)) + c*((1-u).*y1 + u.*((1-v).*y2 + v.*y3))).*exp(-J*k_0*argexp);
f = p.*jacob;

FUNCTION CALL:

fun = @(u,v) abcintegrand(u,v,a(i),b(i),c(i),k_0,x1,x2,x3,y1,y2,y3);
b_integral(e,i) = (E_0*k_0*k_0*(eps_r(e)-1/mu_r(e))/(2*D(e)))*dblquad(fun,0,1,0,1);

"Vivek Saxena" <maverick280857@yahoo.com> wrote in message <ho427a$7d4$1@fred.mathworks.com>...

> ..... when I did the integral, the values all came out wrong .....

  You will find that making a transformation of variables in evaluating multiple integrals using a Jacobian determinant is one of the fundamental methods in advanced calculus. It has been around a long time and you can surely count on it being valid. Note however that you must use the absolute value of the Jacobian determinant! I notice that you did not use the absolute value in your example. I should have made that clear earlier in this thread. See these three links:

 http://en.wikipedia.org/wiki/Integration_by_substitution
 http://wapedia.mobi/en/Integration_by_substitution
 http://planetmath.org/encyclopedia/ChangeOfVariablesInIntegralOnMathbbRn.html

  In your particular case of integration over a triangularly-shaped area, the transformation I gave you has a Jacobian determinant of:

 u*det([x1,x2,x3;y1,y2,y3;1,1,1])

(Try simplifying your expression for 'jacob' and you will see.) However, as I said, what you need is the absolute value of the Jacobian as a factor, so if your three triangle vertices are oriented in clockwise order around the triangle, you would need to change the sign of your 'jacob' variable. You need its absolute value.

  It is obvious that if your original integrand is unity, then its double integral over the triangle will of course be its area. The above Jacobian's absolute value is u multiplied by twice the triangle's area, so the integral in the u,v space will be int, u=1 to u=1, 2*area*u du = area also. The two integrals would thus be equal.

  The same should be true of any integrand you might choose, including the one you used. I would suggest you try experimenting with very simple integrands to begin with, such as f(x,y) = x, y, x*y, or x^2+y^2 - something like that - to build up confidence in the method. With a little effort you should be able to find analytic expressions for both forms of the double integral in these simple cases and be able to show them to be identical.

  If the problem is not due to a negative Jacobian factor, you should realize that dblquad works with a default error tolerance of only six decimal places unless you specify a smaller tolerance. How far off were your results?

> I'm unable to convince myself that the nonlinear transformation you wrote really does map a triangle into a square.

  The inverse of the transformation

 x = (1-u)*x1 + u*((1-v)*x2 + v*x3)
 y = (1-u)*y1 + u*((1-v)*y2 + v*y3)

is this (courtesy of matlab symbolic toolbox)

 t = ((x-x1)*(y3-y2)-(y-y1)*(x3-x2)); % Temporary variable
 u = t/det([x1,x2,x3;y1,y2,y3;1,1,1]);
 v = (y*(x2-x1)-x*(y2-y1)+x1*y2-x2*y1)/t;

  Therefore the following code should demonstrate that this inverse does indeed transform all of the interior of a triangle in x,y into the interior of the unit cube in u,v. I have taken the liberty of using my FEX program, randfixedsum, to obtain barycentric coordinates for creating points within a random triangle. It is available at

 http://www.mathworks.com/matlabcentral/fileexchange/9700

  Here is the code:

 x1 = rand-1.1; y1 = rand; % Choose a random triangle off to the left
 x2 = rand-1.1; y2 = rand;
 x3 = rand-1.1; y3 = rand;
 r = randfixedsum(3,2000,1,0,1); % All 2000 columns sum to 1
 a = r(1,:); b = r(2,:); c = r(3,:);
 x = a*x1+b*x2+c*x3; % a,b,c are barycentric coord.
 y = a*y1+b*y2+c*y3;
 t = ((x-x1)*(y3-y2)-(y-y1)*(x3-x2)); % Here is the inverse transform
 u = t/det([x1,x2,x3;y1,y2,y3;1,1,1]);
 v = (y*(x2-x1)-x*(y2-y1)+x1*y2-x2*y1)./t;
 plot(x,y,'y.',[x1,x2,x3,x1],[y1,y2,y3,y1],'b-',...
      u,v,'r.',[0,1,1,0,0],[0,0,1,1,0],'b-')
 axis equal

  As you can see, it transforms uniformly distributed points within the triangle over to points within the unit square in u,v, though no longer uniformly spaced. This non-uniformity reflects the variation of the Jacobian factor - the density is smaller as u becomes smaller on the left side.

Roger Stafford

Hi Roger, thanks for the detailed reply. I am familiar with the theory of the method (being an engineering student myself :-) ). That is not the problem. I was wondering however, that given the forward transformation (x, y) --> (u, v) which clearly maps the triangle into the square (simply because u and v are now independent), can we show that (u,v) --> (x,y) is unique. Anyway since my last post, I managed to convince myself about this. That was just my way of trying to motivate the choice of the transformation.

Also, when I said that the values of the integral being computed are turning out to be wrong, I was referring not to the transformation but rather to my usage of dblquad. I'm inclined to believe that something is wrong with my usage of that function. I should've been more explicit. (Thanks for your help btw.)

"Vivek Saxena" <maverick280857@yahoo.com> wrote in message <ho6g8s$os1$1@fred.mathworks.com>...
> Hi Roger, thanks for the detailed reply. I am familiar with the theory of the method (being an engineering student myself :-) ). That is not the problem. I was wondering however, that given the forward transformation (x, y) --> (u, v) which clearly maps the triangle into the square (simply because u and v are now independent), can we show that (u,v) --> (x,y) is unique. Anyway since my last post, I managed to convince myself about this. That was just my way of trying to motivate the choice of the transformation.
>
> Also, when I said that the values of the integral being computed are turning out to be wrong, I was referring not to the transformation but rather to my usage of dblquad. I'm inclined to believe that something is wrong with my usage of that function. I should've been more explicit. (Thanks for your help btw.)

One way to check the correctness of the parametrization is to integrate a polynomial P(x) on the triangle. If P(x) the polynomial of
- order 0: the integral result is P(x1)*|T|
- order 1; the integral result is (P(x1)+P(x2)+P(x3)) * |T|/3
etc...

x1, x2, x3 are three corners, and |T| is the area of triangle

Once the implementation gives correct result, you can then replace P(x) by your function.

Bruno

"Vivek Saxena" <maverick280857@yahoo.com> wrote in message <ho6g8s$os1$1@fred.mathworks.com>...
> .........
> Also, when I said that the values of the integral being computed are turning out to be wrong, I was referring not to the transformation but rather to my usage of dblquad. I'm inclined to believe that something is wrong with my usage of that function. .......
------------
  One thing I forgot to ask you, Vivek. Are your values of k_0 and the ranges of x in the triangle such that an excessive number of cycles occur in the trigonometric function, exp(-i*k_0*x)? If so, numerical integration is a correspondingly difficult process and dblquad, or possibly even quad2d, may very well end up confused and erroneous unless you place a very tight tolerance setting in them. This is true whether or not you have performed a transformation of the variables, and it is also true even for single integration.

  Quadrature routines can be temperamental beasts when faced with rapidly fluctuating integrands. You would do well to read the description of quad2d at

 http://www.mathworks.com/access/helpdesk/help/techdoc/ref/quad2d.html

to get an understanding of these difficulties even if you are only using dblquad.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <ho81p7$1nv$1@fred.mathworks.com>...

> One thing I forgot to ask you, Vivek. Are your values of k_0 and the ranges of x in the triangle such that an excessive number of cycles occur in the trigonometric function, exp(-i*k_0*x)?

Roger, I can answer this question. Usually not: when the FEM is correctly setup then the mesh size must be small compare to wave number, thus exp(-i*k_0*x) would not have many cycles on the triangle element. If it's not the case then there is no use to solve FEM because the discretized solution will have huge error.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ho82fo$eg4$1@fred.mathworks.com>...
> "Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <ho81p7$1nv$1@fred.mathworks.com>...
>
> > One thing I forgot to ask you, Vivek. Are your values of k_0 and the ranges of x in the triangle such that an excessive number of cycles occur in the trigonometric function, exp(-i*k_0*x)?
>
> Roger, I can answer this question. Usually not: when the FEM is correctly setup then the mesh size must be small compare to wave number, thus exp(-i*k_0*x) would not have many cycles on the triangle element. If it's not the case then there is no use to solve FEM because the discretized solution will have huge error.
>
> Bruno

Hi Roger and Bruno. Thanks for your replies and help so far. I'm using a unit wavelength, and k_0 = 2*pi as a result. Also, the mesh size has been kept small enough -- actually this is a textbook problem but I have to solve it with a more complicated boundary condition than is given in the standard books. So these particular double integrals arise in the second order absorbing boundary condition, but not in the formulation with the first order absorbing boundary condition (for which the code runs perfectly). In fact, theory tells me that the second order ABC should work fine with a coarser discretization as compared to the first order ABC (if things are going fine). But I am being conservative and keeping the mesh size just as fine as it was earlier.

Roger, I think it is likely that the program is plagued by some numerical errors, as you have pointed out. (I have gone through the code and ensured that there are no logical or syntactic errors.) I must admit that I do not know these techniques very well theoretically, but I will try and play with the options. If you have more suggestions on how the tolerances could be set in view of the kind of function that the integrand is, please do advise me.

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ho75n5$b42$1@fred.mathworks.com>...
> One way to check the correctness of the parametrization is to integrate a polynomial P(x) on the triangle. If P(x) the polynomial of
> - order 0: the integral result is P(x1)*|T|
> - order 1; the integral result is (P(x1)+P(x2)+P(x3)) * |T|/3
> etc...
>
> x1, x2, x3 are three corners, and |T| is the area of triangle
>
> Once the implementation gives correct result, you can then replace P(x) by your function.
>
> Bruno
------------
  Hi Bruno. Your "etc..." intrigued me, so I worked out what the "order 2" test would be. I wasn't aware of this identity before. What would an order 3 test be?

  The double integral of any quadratic polynomial,

 P(x,y) = A*x^2+B*x*y+C*y^2+D*x+E*y+F ,

over a triangle is equal to the triangle's area, multiplied by one quarter of the polynomial's average value at the three vertices plus three quarters of its value at the triangle's centroid (located at the averages of the three respective coordinates.)

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <hoas8r$jbd$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ho75n5$b42$1@fred.mathworks.com>...
> > One way to check the correctness of the parametrization is to integrate a polynomial P(x) on the triangle. If P(x) the polynomial of
> > - order 0: the integral result is P(x1)*|T|
> > - order 1; the integral result is (P(x1)+P(x2)+P(x3)) * |T|/3
> > etc...
> >
> > x1, x2, x3 are three corners, and |T| is the area of triangle
> >
> > Once the implementation gives correct result, you can then replace P(x) by your function.
> >
> > Bruno
> ------------
> Hi Bruno. Your "etc..." intrigued me, so I worked out what the "order 2" test would be. I wasn't aware of this identity before. What would an order 3 test be?
>
> The double integral of any quadratic polynomial,
>
> P(x,y) = A*x^2+B*x*y+C*y^2+D*x+E*y+F ,
>
> over a triangle is equal to the triangle's area, multiplied by one quarter of the polynomial's average value at the three vertices plus three quarters of its value at the triangle's centroid (located at the averages of the three respective coordinates.)
>
> Roger Stafford

Hi Roger,

The order 3 polynomial, the formula is

I = 9/20*T*[P((x1+x2+x3)/3)] +
2/15*T*[P((x1+x2)/2)+P((x2+x3)/2)+P((x3+x1)/2)] +
1/20*T*[P(x1) + P(x2) + P(x3)].

Those kind of Gauss's formula are derived mostly for FEM (I copy it from a book).

Bruno