The following describes and calculates a number of useful properties of an object defined by a triangular surface mesh. Specifically, we consider:
, where is an infinitesimal area on the surface
The volume is given by If the density is uniform, mass is just density times volume.
Assuming its a uniformly dense solid (as opposed to a thin shell, e.g.) the centroid is defined as
, where
In tensor notation we have
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:
where is the mass of the object and is the centroid.
This is defined using tensor notation as
The trace of will be equal to the surface area of the mesh.
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.
, where is the area of the 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.
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.
The surface integral can then be written as a summation of the integral evaluated per 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.
For obtaining the integral that gives us the volume, the function has to be such that , so that the right hand side in the Divergence Theorem equation equals volume. For example, it can be taken to be . In the code I take where , 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
where lies on the surface, and are the components of the outward facing normal at . 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.
x
def 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_tensorif __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)