How to integrate a function over a 3D planar polygon?

I am trying to integrate a function over a planar polygon in 3D. In 2D, this is fairly straightforward, using either answer from this question (I use the second answer). If we use an equilateral triangle centered at the origin as our polygon:

f[pos_List] := 1
basetri={{1/2,- (1/6) Sqrt[3],0},{0,1/3 Sqrt[3],0},{-(1/2),- (1/6) Sqrt[3],0}}
inPolyQ[poly_,pt_]:=2.π==Total[VectorAngle@@@Transpose@{#,RotateRight[#]}]&@(#-pt&/@poly)


we can integrate our function f (in this case, this function will return the area of the triangle):

NIntegrate[
f[{x,y,0}] Boole[inPolyQ[basetri,{x,y,0}]],
{x,Min@basetri[[;;,1]],Max@basetri[[;;,1]]},
{y,Min@basetri[[;;,2]],Max@basetri[[;;,2]]}
]
(* Out[]:= 0.43265 *)


which is pretty close to the exact value of Sqrt[3]/4 (~0.433013). (Yes, the are a bunch of errors thrown.) This can also be done over infinite bounds (see this answer):

NIntegrate[f[{x,y,0}] Boole[inPolyQ[basetri,{x,y,0}]],
{x,-Infinity,Infinity},{y,-Infinity,Infinity}
]
(* Out[]:= 0.433321 *)


which is also pretty close to the exact value. (I'll leave alone for now the question of which bounds should give better results...)

The problem arises when doing this in 3D:

NIntegrate[
f[{x,y,z}] Boole[inPolyQ[basetri, {x, y, z}]],
{x, -Infinity,Infinity}, {y, -Infinity, Infinity}, {z, -Infinity, Infinity}
]
(* Out[]:= 0 *)


It gives this error:

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option.

Increasing MinRecursion doesn't help. I suspect that because this is a planar polygon, it has trouble figuring out exactly where that infinitesimally-thin plane is that satisfies the inequality actually is.

So, how can I integrate a function over an arbitrarily-oriented planar polygon in 3D? For what it's worth, at the very least it would be an acceptable answer for my current problem if it just works with a triangular region (in 3D), but I figured I'd leave the question more general if there's a way to do it.

-

Even if it's a 2D planar polygon in 3D space, it's still a 2D object. 2D objects are always 0 under a 3D measure.

Here we go through two methods of performing the scalar surface integral with a random polygon example.

Say we have a general 3D polygon defined as this:

polygonPts3D = {
{-0.902757, -0.116805, 0},
{0.203504, -0.972294, 0},
{0.849893, 0.414192, 0},
{0.374057, 0.835407, 0},
{-0.907079, 0.352119, 0}
};


and a scalar function defined in $xyz$ space:

Clear[f]
f[x_,y_,z_] := x^3-y^2+z


The standard way:

It's easy to derive the equation of the planar (planeEq) which the polygon is on:

Clear[normFunc]
normFunc[pts_] := Mean[
Normalize[Cross[##]] & @@@
Partition[
Subtract @@@ Partition[pts, 2, 1, 1],
2, 1, 1]
]

norm = normFunc[polygonPts3D]

planeEq = norm.({x, y, z} - Mean[polygonPts3D]) == 0

ContourPlot3D[
Evaluate[planeEq],
{x, -1, 1}, {y, -1, 1}, {z, -1, 1},
RegionFunction -> Function[{x, y, z, func},
inPolyQ[ReplacePart[polygonPts3D, {_, 3} :> 0], {x, y, 0}]
],
Mesh -> 20,
MeshFunctions -> Function[{x, y, z, func}, f[x, y, z]],
ColorFunction -> Function[{x, y, z, func}, Hue[f[x, y, z]]],
ColorFunctionScaling -> False,
AxesLabel -> (Style[#, Red, Bold, 20] & /@ {x, y, z})
]


Thus the standard way to perform the surface integral

$$\int_P f(x,y,z)\,\mathrm{d}S = \int_{P_{xy}} f(x,y,z(x,y)) \sqrt{1+\left(\frac{\partial z}{\partial x}\right)^2+\left(\frac{\partial z}{\partial y}\right)^2} \,\mathrm{d}x\mathrm{d}y$$

would be straightforward in Mathematica

zFunc[x_, y_] := Evaluate[z /. Solve[planeEq, z][[1]]]

1/norm[[3]] NIntegrate[
f[x, y, zFunc[x, y]]*
Boole[inPolyQ[ polygonPts3D, {x, y, zFunc[x, y]} ]],
{x, -∞, ∞}, {y, -∞, ∞}
] // Quiet


0.484606

The other way:

Also, we can establish a new coordinate system $x'y'z'$ such that the polygon lies on plane of $z'=0$, so we can convert the surface integral in $xyz$-space to a 2D area integral in $x'y'$-space.

To achieve that, we need a rotation which convert the norm of the polygon to $z'$ axis:

transMatrix = RotationMatrix[{
normFunc[polygonPts3D],
{0, 0, 1}
}]\[Transpose]
transMatrixInverse = Inverse@transMatrix


So the polygon in $x'y'z'$ should be:

polygonPts2D = ReplacePart[polygonPts3D.transMatrix, {_, 3} :> 0]


and the distance from the origin to the polygon plane:

distance = norm.Mean[polygonPts3D]


It can be seen that the transformed polygon is essentially the same with the original one:

RegionPlot[
Evaluate[inPolyQ[polygonPts2D, {x, y, 0}]],
{x, -1.2, 1.2}, {y, -1.2, 1.2},
MeshFunctions -> Function[{x, y, z, func},
f @@ ({{x, y, distance}}.transMatrixInverse)[[1]]
],
Mesh -> 20,
ColorFunction -> Function[{x, y, z, func},
Hue[f @@ ({{x, y, distance}}.transMatrixInverse)[[1]]]
],
ColorFunctionScaling -> False,
FrameLabel -> (Style[#, Red, Bold, 20] & /@ {x, y})
]


And the integral is:

NIntegrate[
f @@ ({{x, y, distance}}.transMatrixInverse)[[1]]*
Boole[inPolyQ[polygonPts2D, {x, y, 0} ]],
{x, -∞, ∞}, {y, -∞, ∞}
] // Quiet


0.483194

Side note:

For planar polygons, I think there is conformal transformation which can provide bijective maps between the polygon and the unit disk (or eventually maps between the polygon and a rectangle), which in some cases might simplify the integrals and provide more precise results.

-
This works fine for the area function; I guess you'd need to map the {x,y,z} coordinates passed into the function, as well? –  Eli Lansey Jun 21 '13 at 16:03
@EliLansey I think once you have the coordinates of the 3D polygon, finding its norm thus performing the coordinate transformation are both automatically done. –  Silvia Jun 21 '13 at 16:10
I might be misunderstanding something in your answer. If the function is, e.g. f[x_,y_z_]:=x+y+z and your triangle's normal is +y, the function in the integrand needs to be evaluated at the actual x, y z. In my example, f[x_,y_z_]:=1 so there's no need to worry about that, but for an arbitrary function you'd need to map those coordinates, too. –  Eli Lansey Jun 21 '13 at 16:51
@EliLansey OK I get your point. It seems I overlooked your title... But please allow me update the answer after 20 hours. It's too late here and I get some stuff to busy with tomorrow. (And yes, the only thing you need to do is doing the some coordinate transformation on arguments of f.) –  Silvia Jun 21 '13 at 17:38
I've updated my question to make the fact that I'm integrating a function a little clearer. And no rush -- for the time being, my answer below is sufficient for my needs, but not general. –  Eli Lansey Jun 21 '13 at 18:13

For a triangle, this can be solved by integrating in barycentric coordinates. Assuming the function that is being integrated takes in f[{x,y,z}] as an argument:

f[pos_List] := 1
integralOverTriangle[vertices_,f_]:=
Module[{area},
area=Norm[1/2 Cross[vertices[[1]]-vertices[[2]],vertices[[1]]-vertices[[3]]]];
2 area Integrate[f[l1 vertices[[1]]+l2 vertices[[2]]+(1-l1-l2)vertices[[3]]],
{l2,0,1},  {l1,0,1-l2}]
]
nintegralOverTriangle[vertices_,f_]:=
Module[{area},
area=Norm[1/2 Cross[vertices[[1]]-vertices[[2]],vertices[[1]]-vertices[[3]]]];
2 area NIntegrate[f[l1 vertices[[1]]+l2 vertices[[2]]+(1-l1-l2)vertices[[3]]],
{l2,0,1},  {l1,0,1-l2}]
]


Then,

integralOverTriangle[basetri,1&]
(* Out[]:= Sqrt[3]/4 *)


and, the triangle can be oriented arbitrarily:

integralOverTriangle[RotationTransform[{{0, 0, 1}, Normalize[{1, 1, 1}]}] /@ basetri, 1 &]
(* Out[]:= Sqrt[3]/4 *)

-

In Mathematica 10, this can be easily computed using built-in functionality. Using Silvia's answer's polygon and function:

polygonPts3D = {
{-0.902757, -0.116805, 0},
{0.203504, -0.972294, 0},
{0.849893, 0.414192, 0},
{0.374057, 0.835407, 0},
{-0.907079, 0.352119, 0}
};
f[x_,y_,z_] := x^3-y^2+z


All you need to do is:

Integrate[f[x, y, z], Element[{x, y, z}, Polygon[polygonPts3D]]]
Out[]:= -0.323842

-