'------------------------------------------------------------------------------
' NAME:		FGetVolume. Rejean Gagne
'
' OUTPUT:   Volume of the polygon mesh
'
' DESCRIPTION: Calculates the volume (for polygon meshes only).
'
' Method used comes from Graphics Gem V, p.37
'
' Volume = 1/6 * SUM(j=0, m-1) Pface(j) dot (2*Area(j))
' Area = SUM(i=1,h-1) (P2i-P0) X (P2i+1-P2i-1)  + (P2h-P0) X (Pl-P2h-1)
' where
'    Volume (of the mesh) : m = number of facets
'							Pface(j) is an arbitrary vertex of the face (we take the 1st one)
'    Area (of a facet) : k = number of vertices, 
'						 h=1/2*(k-1),
'						 l = 0 if k is odd, l = k-1 if k is even
'------------------------------------------------------------------------------
Set oSel  = Selection(0)
FGetVolume oSel

function FGetVolume(l_meshobj)

	FGetVolume = 0

	l_vol = 0
	l_area = 0

	' We will access to facets as well as the points of the mesh
	set oGeometry = l_meshobj.obj
	set oGeom2D = oGeometry.Geometry2D
	set oGeom0D = oGeometry.Geometry0D

	' We need a bunch of vectors, named after the equations above	
	set oPos  	= CreateObject("Sumatra\Scripting\Math\SIVector3")
	set oP2i  	= CreateObject("Sumatra\Scripting\Math\SIVector3")
	set oP0   	= CreateObject("Sumatra\Scripting\Math\SIVector3")
	set oP2ip1  = CreateObject("Sumatra\Scripting\Math\SIVector3")
	set oP2im1  = CreateObject("Sumatra\Scripting\Math\SIVector3")
	set oPa	= CreateObject("Sumatra\Scripting\Math\SIVector3")
	set oPb	= CreateObject("Sumatra\Scripting\Math\SIVector3")
	set oPc	= CreateObject("Sumatra\Scripting\Math\SIVector3")
	set oParea	= CreateObject("Sumatra\Scripting\Math\SIVector3")
	set oPface	= CreateObject("Sumatra\Scripting\Math\SIVector3")	
	
							
	' For each facet/polygon...
	for i = 0 to oGeometry.Nb2D - 1

		' Get the facet
		o0DSubComp = oGeom2D.SubComponent0D( i )
		
		' intermediate indices
		k = UBound( o0DSubComp, 1 ) + 1
		h = ( k - 1 ) \ 2  ' note the integer division '\'
		l = 0
		if ( k mod 2 ) = 0 then
			l = k - 1
		end if	
		
		' Init the area vector
		oParea.Set 0.0, 0.0, 0.0
		
		'logmessage " Face " & i & " k,h,l = [" & k & "," & h & "," & l & "]" 	

		' position of first vertex is needed
		oGeom0D.Position o0DSubComp( 0 ), oPos  	
		oP0.Set oPos.x, oPos.y, oPos.z

		' this ends up being used for facets with more than 4 vertices		
		for j = 1 to h - 1
		

			oGeom0D.Position o0DSubComp( 2*j ), oPos  	
			oP2i.Set oPos.x, oPos.y, oPos.z
			
			oGeom0D.Position o0DSubComp( 2*j + 1 ), oPos  	
			oP2ip1.Set oPos.x, oPos.y, oPos.z
			
			oGeom0D.Position o0DSubComp( 2*j - 1 ), oPos  	
			oP2im1.Set oPos.x, oPos.y, oPos.z
	
			oPa.Sub oP2i, oP0
			oPb.Sub oP2ip1, oP2im1		
			oPc.Cross oPa, oPb
			
			' cumulate the area	vector		
			oParea.AddInPlace  oPc

			l_area = l_area + oParea.Length
			
		next
		
		' Last component (see equations)
		oGeom0D.Position o0DSubComp( 2 * h ), oPos  	
		oP2i.Set oPos.x, oPos.y, oPos.z
			
		oGeom0D.Position o0DSubComp( l ), oPos  	
		oP2ip1.Set oPos.x, oPos.y, oPos.z
			
		oGeom0D.Position o0DSubComp( 2*h - 1 ), oPos  	
		oP2im1.Set oPos.x, oPos.y, oPos.z
	
		oPa.Sub oP2i, oP0
		oPb.Sub oP2ip1, oP2im1		
		oPc.Cross oPa, oPb
		
		' cumulate the area	vector			
		oParea.AddInPlace  oPc
		
		' while we're at it, cumulate the surface area		
		l_area = l_area + oParea.Length

		' compute the voume contribution for this facet
		oGeom0D.Position o0DSubComp( 0 ), oPface
	
		l_dot = oPface.Dot( oParea )				
		l_volume = l_volume + l_dot
	
	next
	
	' Final volume
	FGetVolume = l_volume / 6

	logmessage "volume " & FGetVolume

end function