# Calculating properties of solids represented by triangular surface mesh

The following describes and calculates a number of useful properties of an object defined by a triangular surface mesh. Specifically, we consider:

• Surface Area
• Volume ($\propto$ mass, $M$, assuming uniform density)
• Centroid
• Moment of Inertia
• Area tensor or Surface energy tensor

## Definitions

#### Surface Area

$A_{total} = \int dA$, where $dA$ is an infinitesimal area on the surface

#### Volume

The volume is given by ​ $V = \int dV$ If the density is uniform, mass is just density times volume.

#### Centroid

Assuming its a uniformly dense solid (as opposed to a thin shell, e.g.) the centroid is defined as $\bar{r}_c = \dfrac{\int \rho~\bar{x} dV}{\int \rho~dV} = \dfrac{\int \bar{x} dV}{\int dV}$

#### Moment of Inertia/Inertia Tensor

$\mathbf{I} = \begin{bmatrix} I_{xx} & I_{xy} & I_{xz} \\ I_{xy} & I_{yy} & I_{yz} \\ I_{xz} & I_{yz} & I_{zz} \end{bmatrix}$, where

$I_{xx} = \int\rho (y^2+z^2)dV$

$I_{yy} = \int\rho (x^2+z^2)dV$

$I_{zz} = \int\rho (x^2+y^2)dV$

$I_{xy} = -\int\rho xy~dV$

$I_{xz} = -\int\rho xz~dV$

$I_{yz} = -\int\rho yz~dV$

In tensor notation we have ​ $\mathbf{I}_{ij} = \int{\rho (x_k x_k\delta_{ij} - x_i x_j)dV }$

In general the above integral give the moment of inertia about axes passing through the origin of the coordinate system. To get the body centered moment of inertia we can use the parallel-axis theorem: $\mathbf{I}^{BC} = \begin{bmatrix} I_{xx} & I_{xy} & I_{xz} \\ I_{xy} & I_{yy} & I_{yz} \\ I_{xz} & I_{yz} & I_{zz} \end{bmatrix} - M\begin{bmatrix} y_c^2+z_c^2 & -x_c y_c & -x_c z_c \\ -x_c y_c & x_c^2+z_c^2 & -y_c z_c \\ -x_c z_c & -y_c z_c & x_c^2+y_c^2 \end{bmatrix}$

where $M$ is the mass of the object and $(x_c,y_c,z_c)$ is the centroid.

### Area Tensor/Surface Energy Tensor

This is defined using tensor notation as

$\mathbf{S}_{ij} = \dfrac{1}{2}\int(\delta_{ij}- n_i n_j)dA$

The trace of $\mathbf{S}$ will be equal to the surface area of the mesh.

## Computational Details

This only works if all the triangle facets have to be oriented consistently. That is, in the face data, for each face the vertices have to be listed in a counter-clockwise order, while looking at the face from outside the mesh. If the vertices are consistently clockwise, that would result in a negative value for volume. In such a case, switching any two columns in the face data will make it oriented properly. The surface integrals are evaluated exactly directly on the triangles. So, the total surface area and the area tensor can be calculated in each face and added up.

### Surface Area

$A_{total} = \int dA = \sum\limits_{i}^{}{A_{i}}$, where $A_i$ is the area of the $i^\text{th}$ face. The area of each triangle is calculated by taking cross product of two edges of a triangle. Half the magnitude of this vector is the area of the triangle. This step also gives us the normal to the triangle face, which will be needed later on.

### Volume integrals

The volume integrals can be done directly by decomposing the mesh into tetrahedrons and evaluating the integral for each tetrahedron. An elegant alternative is to use the Gauss Divergence Theorem. We can use it to recast the volume integrals as surface integrals.

$\int(\nabla\cdot\bar f)dV = \int(\bar f\cdot \bar n)dA$

The surface integral can then be written as a summation of the integral evaluated per face.

$\int(\bar f\cdot \bar n)dA =\sum\limits_{\text{faces}} \int(\bar f\cdot \bar n)dA_{\textit{face}}$

The integrand is calculated exactly on each triangular face. There are 10 functions whose volume integrals need to be evaluated-- 1 for volume, 3 for centroid and 6 for the moment of inertia.

#### Volume

For obtaining the integral that gives us the volume, the function $f$ has to be such that $\nabla\cdot\bar f=1$, so that the right hand side in the Divergence Theorem equation equals volume. For example, it can be taken to be $x \mathbf{e_x}+0\mathbf{e_y}+0\mathbf{e_z}$. In the code I take $f(x)=\dfrac{1}{3}\bar{r}$ where $\bar{r} = (x\mathbf{e_x}+y\mathbf{e_y}+z\mathbf{e_z})$, whose divergence is also unity. This is done because this particular function is needed later in calculating the other integrals.

Plugging the latter into the equation of the divergence theorem, we get

$\int dV =\int(\nabla\cdot(\dfrac{1}{3}(x\mathbf{e_x}+y\mathbf{e_y}+z\mathbf{e_z})))dV$

$~~~~~~~~= \int(\dfrac{1}{3}(x\mathbf{e_x}+y\mathbf{e_y}+z\mathbf{e_z}))\cdot \bar n)dA$

$~~~~~~~~=\dfrac{1}{3}\int(\bar{r}\cdot\bar{n})dA = \dfrac{1}{3}\int(xn_x+yn_y+zn_z)dA$

where $\bar{r}$ lies on the surface, and $(n_x,n_y,n_z)$ are the components of the outward facing normal at $\bar{r}$. As mentioned earlier, the surface integral can be written as sum of exact per-face integrals.

The others are calculated similarly. The python code follows. The sample program calculates the the mass properties for triangulated unit cube of unit density.

xdef get_mesh_props(vertices,faces):    '''    Abhilash Reddy Malipeddi. abhilashreddy.com    Calculates volume, centroid and centered second moment of inertia tensor    for a trimesh. The corner points of each triangle have to be numbered     counter-clockwise(CCW), while looking from outside the mesh. The volume    calculated by this code will be positive if the vertices are oriented correctly.     '''    import numpy as np​    # Pick new origin close to the action, i.e. near the mesh to reduce     # roundoff errors in case the mesh is located far away from the origin.    # Doing this the following way to keep the original array untouched    vertices1=vertices.copy()    ref=np.mean(vertices1,axis=0) # temporary origin for calculation    vertices1-=ref ​    # Face centroids for each triangle      cnt=(1.0/3.0)*(vertices1[faces[:,0]]+vertices1[faces[:,1]]+vertices1[faces[:,2]])​    # triangle edge vectors    e0=vertices1[faces[:,1]]-vertices1[faces[:,0]]    e1=vertices1[faces[:,2]]-vertices1[faces[:,0]]​    FN=0.5*np.cross(e0,e1)                         # triangle/face normal vector with magnitude equal to area of the face     ar=np.sqrt(np.sum(FN**2,axis=1))               # Surface area of each triangle    FUN=FN/ar[:,np.newaxis]                        # Face unit outward normal    vol  = np.sum(cnt*FUN*ar[:,np.newaxis])/3.0    # Volume of the mesh    area = np.sum(ar)                              # Surface area of mesh​    C2F=cnt**2*FUN*ar[:,np.newaxis]​    centroid=np.sum(C2F,axis=0)*0.5/vol    xc,yc,zc=centroid​    pxx = - vol*xc**2 + 1.0/3.0*(np.sum(cnt[:,0]*C2F[:,0]))     pyy = - vol*yc**2 + 1.0/3.0*(np.sum(cnt[:,1]*C2F[:,1]))     pzz = - vol*zc**2 + 1.0/3.0*(np.sum(cnt[:,2]*C2F[:,2]))     pxy =   vol*xc*yc - 1.0/4.0*(np.sum(cnt[:,0]*C2F[:,1] + cnt[:,1]*C2F[:,0]))     pyz =   vol*yc*zc - 1.0/4.0*(np.sum(cnt[:,1]*C2F[:,2] + cnt[:,2]*C2F[:,1]))     pzx =   vol*zc*xc - 1.0/4.0*(np.sum(cnt[:,2]*C2F[:,0] + cnt[:,0]*C2F[:,2])) ​    inertia_tensor = np.reshape([ pyy+pzz,      pxy,      pzx,                                      pxy,  pxx+pzz,      pyz,                                      pzx,      pyz,  pxx+pyy  ],(3,3))​    Axx =  np.sum( (1.0 - FUN[:,0]**2      )*ar)    Ayy =  np.sum( (1.0 - FUN[:,1]**2      )*ar)    Azz =  np.sum( (1.0 - FUN[:,2]**2      )*ar)    Axy =  np.sum( (    - FUN[:,0]*FUN[:,1])*ar)    Axz =  np.sum( (    - FUN[:,0]*FUN[:,2])*ar)    Ayz =  np.sum( (    - FUN[:,1]*FUN[:,2])*ar)​    area_tensor = 0.5*np.reshape([ Axx, Axy, Axz ,                               Axy, Ayy, Ayz ,                               Axz, Ayz, Azz  ],(3,3))​# add the offest back to the centroid    centroid+=ref    return area,vol,centroid,area_tensor,inertia_tensor​if __name__ == "__main__":      import numpy as np# Output:#    Area:  6.0 , Volume:  1.0 , Centroid:  [0.5 0.5 0.5]#    Moment of Inertia:#     [[0.16666667 0.         0.        ]#     [0.         0.16666667 0.        ]#     [0.         0.         0.16666667]]#    Area Tensor:#     [[2. 0. 0.]#     [0. 2. 0.]#     [0. 0. 2.]]      vertices=np.asarray([               [0,0,0],               [0,0,1],               [0,1,0],               [1,0,0],               [0,1,1],               [1,0,1],               [1,1,0],               [1,1,1],],dtype=np.float64)​      faces=np.asarray([                  [2,4,7],                  [2,7,6],                  [2,6,3],                  [2,3,0],                  [7,5,3],                  [7,3,6],                  [4,1,7],                  [7,1,5],                  [0,5,1],                  [0,3,5],                  [2,0,4],                  [4,0,1],],dtype=np.int)​      area,volume,centroid,area_tensor,I = get_mesh_props(vertices,faces)      print('Area: ',area,', Volume: ',volume,', Centroid: ',centroid)      print('Moment of Inertia: \n',I)      print('Area Tensor: \n', area_tensor)