Lately I was thinking about computing section properties programatically. Of course, any engineering undergrad can compute the area and moment of inertia of a section comprising rectangles and other simple shapes — for example, a tee-section or an I-beam section. There is software out there that computes section properties of arbitrary sections, including the area, centroid and moment of inertia. Presumably that software is doing something other than just breaking down the shape into a bunch of rectangles or other simple shapes.

As I was thinking about this, I remembered back nearly twenty years to an undergrad course that taught Green’s Theorem. This is a way to compute an area integral using a contour intergral. It seems to me that computing a contour integram of an arbitrary shape is probably easier than computing an area integral.

Green’s theorem is as follows:

$$\oint_C \left(L\,dx + M\,dy\right) = \iint_D \left( \frac{\partial M}{\partial x} - \frac{\partial L}{\partial y} \right)\,dA$$

Area

The first section property that is obviously needed is area. Area is (clearly) given by:

$$Area = \iint_D 1\,dA$$

This means that to apply Green’s theorum, we need to select functions \(M\) and \(L\) such that:

$$\frac{\partial M}{\partial x} - \frac{\partial L}{\partial y} = 1$$

There are an infinite number of functions that would satisfy that equality, but let’s choose:

$$M = \frac{1}{2} x \\ L = -\frac{1}{2} y$$

We’ll need to integrate these functions over the closed contour enclosing the shape. It will be convinient to evaluate that contour interval if the shape can be described as piecewise functions of the parameter \(s\) that varies from \(0\) to \(1\) over each piecewise segment describing the closed contour. With this in mind, we’ll define two functions as follows:

$$f\left(s\right) = x \\ g\left(s\right) = y$$

For a linear segment with endpoints \(\left(x_a, y_a\right)\) and \(\left(x_b, y_b\right)\), these functions would be as follows. Similar functions can be written for other segments.

$$x = f\left(s\right) = x_a + \left(x_b - x_a\right) s \\ y = g\left(s\right) = y_a + \left(y_b - y_a\right) s$$

Noting the following:

$$f^\prime\left(s\right)\,ds = dx \\ g^\prime\left(s\right)\,ds = dy$$

We can make the following substitutions:

$$Area = \iint_D 1\,dA \\ = \oint_C \left(L\,dx + M\,dy\right) \\ = \oint_C \left(\frac{-1}{2} y\,dx + \frac{1}{2} x\,dy\right) \\ = \sum \int \left(\frac{-1}{2} g_i\left(s\right)\,dx + \frac{1}{2} f_i\left(s\right)\,dy\right) \\ = \sum \int \left(\frac{-1}{2} g_i\left(s\right)\,f_i^\prime\left(s\right)\,ds + \frac{1}{2} f_i\left(s\right)\,g_i^\prime\left(s\right)\,ds\right) \\ = \frac{1}{2} \sum\int \left[f_i\left(s\right)\,g_i^\prime\left(s\right) - g_i\left(s\right)\,f_i^\prime\left(s\right)\right]\,ds$$

Example Implementation of Area Computation

For a simple example, let’s consider the computation of the area bounded by a set of line segments. The functions \(f\) and \(g\) for this case are given above; the derivatives, which are also needed, are as follows:

$$f^\prime\left(s\right) = \left(x_b - x_a\right) \\ g^\prime\left(s\right) = \left(y_b - y_a\right)$$

The area bounded by \(n\) such line segments is:

$$Area = \frac{1}{2} \sum_{i=1}^{n} \int \left[f_i\left(s\right)\,g_i^\prime\left(s\right) - g_i\left(s\right)\,f_i^\prime\left(s\right)\right]\,ds \\ = \frac{1}{2} \sum_{i=1}^{n} \int \left\{ \left[x_{ai} + \left(x_{bi} - x_{ai}\right) s\right]\left[y_{bi} - y_{ai}\right] - \left[y_{ai} + \left(y_{bi} - y_{ai}\right) s\right]\left[x_{bi} - x_{ai}\right] \right\}\,ds \\ = \frac{1}{2} \sum_{i=1}^{n} \int \left\{ x_{ai}\left(y_{bi} - y_{ai}\right) - y_{ai}\left(x_{bi} - x_{ai}\right) \right\}\,ds \\ = \frac{1}{2} \sum_{i=1}^{n} \left\{ x_{ai}\left(y_{bi} - y_{ai}\right) - y_{ai}\left(x_{bi} - x_{ai}\right) \right\}$$

This can be implemented in Python. This implementation will take a numpy array of size (N, 2) representing the vertecies.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches

def compute_area(v):
  v_offset = np.concat((v[1:], v[0:1]))
  return 1 / 2 * float(sum([
    xa * (yb - ya) - ya * (xb - xa)
    for (xa, ya), (xb, yb) in zip(v, v_offset)
  ]))

We’ll also create a function that plots the shape.

def plot_shape(v):
  fig, ax = plt.subplots()
  ax.add_patch(matplotlib.patches.Polygon(v))
  ax.autoscale(axis="both")

Let’s try a few examples. First, a 4x2 rectangle centered on the origin.

rect2_4 = np.array([
  [-2, -1],
  [2, -1],
  [2, 1],
  [-2, 1]
])

plot_shape(rect2_4)

The area should be \(4 \times 2 = 8\). Let’s check what the function returns:

compute_area(rect2_4)

8.0

Next, let’s try a triangle.

tri3_2 = np.array([
  [0, 0],
  [3, 0],
  [0, 2]
])

plot_shape(tri3_2)

The area should be \(\frac{1}{2} \times 3 \times 2 = 3\). Let’s check what the function returns:

compute_area(tri3_2)

3.0

Computation of the Centroid

The centroid in the \(X\) axis is given by:

$$\bar{X} = \frac{1}{A} \iint_D x\,dA$$

Using Green’s theorem again, we need to choose two new functions \(M\) and \(L\) such that

$$\frac{\partial M}{\partial x} - \frac{\partial L}{\partial y} = x$$

We’ll choose:

$$M = \frac{1}{2} x^2 \\ L = 0$$

As before, we’ll use the functions \(f\left(s\right)\) and \(g\left(s\right)\)

$$\bar{X} = \frac{1}{A} \iint_D x\,dA \\ = \frac{1}{A} \oint_C \left(L\,dx + M\,dy\right) \\ = \frac{1}{A} \oint_C \left(0\,dx + \frac{1}{2} x^2\,dy\right) \\ = \frac{1}{2A} \sum\int \left[f_i\left(s\right)\right]^2\,g_i^\prime\left(s\right)\,ds$$

When the shape comprises \(n\) line segments, this becomes:

$$\bar{X} = \frac{1}{2A} \sum\int \left[x_{ai} + \left(x_{bi} - x_{ai}\right) s\right]^2\,\left(y_{bi} - y_{ai}\right)\,ds \\ = \frac{1}{2A} \sum_int \left[x_{ai}^2 + 2 x_{ai} \left(x_{bi} - x_{ai}\right) s + \left(x_{bi} - x_{ai}\right)^2 s^2\right]\,\left(y_{bi} - y_{ai}\right)\,ds \\ = \frac{1}{2A} \sum_{i=1}^n \left[x_{ai}^2 + x_{ai} \left(x_{bi} - x_{ai}\right) + \frac{1}{3}\left(x_{bi} - x_{ai}\right)^2\right]\,\left(y_{bi} - y_{ai}\right)$$

The centroid in the \(Y\) axis can be found the same way, but reversing the \(x\)‘s and \(y\)‘s and reversing the sign.

Example Implementation of Centroid Computation

This can be implemented in Python as follows:

def compute_centroid(v):
  area = compute_area(v)
  v_offset = np.concat((v[1:], v[0:1]))
  x_bar =  1 / (2 * area) * float(sum([
    (xa**2 + xa * (xb - xa) + 1 / 3 * (xb - xa)**2) * (yb - ya)
    for (xa, ya), (xb, yb) in zip(v, v_offset)
  ]))
  y_bar =  -1 / (2 * area) * float(sum([
    (ya**2 + ya * (yb - ya) + 1 / 3 * (yb - ya)**2) * (xb - xa)
    for (xa, ya), (xb, yb) in zip(v, v_offset)
  ]))
  return x_bar, y_bar

Trying this with the previous example of a 4x2 rectangle centered on the origin should result with the centroid on the origin.

compute_centroid(rect2_4)

(0.0, -0.0)

For the previous triangle, the centroid should be located \(1/3\) of the way from the base.

compute_centroid(tri3_2)

(1.0, 0.6666666666666666)

Computation of the Moment of Inertia

When computing the moment of inertia, the origin needs to be set coincident with the centroid. In the following derrivation, it will be assumed that the centroid is coincident with \(\left(x, y\right) = \left(0, 0\right)\). The moment of inertia about the \(Y\) asis is:

$$I_{YY} = \iint_D x^2\,dA$$

Using Green’s theorem again, we need to choose two new functions \(M\) and \(L\) such that

$$\frac{\partial M}{\partial x} - \frac{\partial L}{\partial y} = x^2$$

We’ll choose:

$$M = \frac{1}{3} x^3 \\ L = 0$$

As before, we’ll use the functions \(f\left(s\right)\) and \(g\left(s\right)\)

$$I_{YY} = \iint_D x^2\,dA \\ = \oint_C \left(L\,dx + M\,dy\right) \\ = \oint_C \left(0\,dx + \frac{1}{3} x^3\,dy\right) \\ = \frac{1}{3} \sum\int \left[f_i\left(s\right)\right]^3\,g_i^\prime\left(s\right)\,ds$$

When the shape comprises \(n\) line segments, this becomes:

$$I_{YY} = \frac{1}{3} \sum\int \left[x_{ai} + \left(x_{bi} - x_{ai}\right) s\right]^3\,\left(y_{bi} - y_{ai}\right)\,ds \\ = \frac{1}{3} \sum\int \left[x_{ai}^3 + 3 x_{ai}^2 \left(x_{bi} - x_{ai}\right) s + 3 x_{ai} \left(x_{bi} - x_{ai}\right)^2 s^2 + \left(x_{bi} - x_{ai}\right)^3 s^3\right]\,\left(y_{bi} - y_{ai}\right)\,ds \\ = \frac{1}{3} \sum_{i=1}^n \left[x_{ai}^3 + \frac{3}{2} x_{ai}^2 \left(x_{bi} - x_{ai}\right) + x_{ai} \left(x_{bi} - x_{ai}\right)^2 + \frac{1}{4}\left(x_{bi} - x_{ai}\right)^3\right]\,\left(y_{bi} - y_{ai}\right)$$

The moment of inertia in the \(X\) axis can be found the same way, but reversing the \(x\)‘s and \(y\)‘s and reversing the sign.

Example Implementation of Moment of Inertia Computation

This can be implemented in Python as follows:

def compute_moi(v):
  centroid = compute_centroid(v)
  v = v - centroid
  v_offset = np.concat((v[1:], v[0:1]))
  Iyy =  1 / 3 * float(sum([
    (xa**3 + 3 / 2 * xa**2 * (xb - xa) + xa * (xb - xa)**2 + 1 / 4 * (xb - xa)**3) * (yb - ya)
    for (xa, ya), (xb, yb) in zip(v, v_offset)
  ]))
  Ixx =  -1 / 3 * float(sum([
    (ya**3 + 3 / 2 * ya**2 * (yb - ya) + ya * (yb - ya)**2 + 1 / 4 * (yb - ya)**3) * (xb - xa)
    for (xa, ya), (xb, yb) in zip(v, v_offset)
  ]))
  return Ixx, Iyy

Trying this with the previous example of a 4x2 rectangle. The moments of inertia should be:

$$I_{xx} = \frac{1}{12} \left(4\right) \left(2\right)^3 = 2.667 \\ I_{yy} = \frac{1}{12} \left(2\right) \left(4\right)^3 = 10.667$$
compute_moi(rect2_4)

(2.6666666666666665, 10.666666666666666)

For the previous 3x2 right triangle, the moment of inertia should be:

$$I_{xx} = \frac{1}{36} \left(3\right) \left(2\right)^3 = 0.667 \\ I_{yy} = \frac{1}{36} \left(2\right) \left(3\right)^3 = 1.500$$
compute_moi(tri3_2)

(0.6666666666666666, 1.5)

Conclusion

Using Green’s Theorem can be used to help compute section properties of shapes given the contour of the shape. In this post, we computed the area, centroid location and moment of inertia of arbitrary closed shapes bounded by a series of line segments. This approach could be extended to conisder other bounding curves (for example, circular arcs, etc.). While I’m not sure how commercial software computes section properties of arbitrary sections, this approach seems like a pretty good guess.