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

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.

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

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