Bending Moment Extraction

The following example demonstrates how internal stresses computed using Scan&Solve can be integrated to approximate the bending moment acting in a beam model. A short video showing the script in operation may be viewed here.

function to numerically integrate SnS solution components (lines 13-42)
function to compute triangle area (lines 46-51)
function to integrate stress over a trimmed surface (lines 57-90)
the main routine (lines 93-130)

Option Explicit
'Script written by Michael Freytag
'Script version Monday, September 16, 2013 3:10:56 PM

'variable for the SnSScript plug-in object
Public objSnSPlugin

'quadIntegrateMoments: Integrate over the triangle defined by the array of vertices 
'in arrVerts using a 7 point quadrature rule
'strObject -the solid with SnS solution data
'arrVerts -the vertices of the triangle to integrate over
'arrCentroid -the centroid of the surface around which moments are being computed
Function IntegrateMoments(strObject, arrVerts, arrCentroid)
        Dim quadWts:quadWts = Array(0.225, 0.13239415, 0.13239415, 0.13239415, 0.12593918, 0.12593918, 0.12593918)
        Dim quadCoords:quadCoords = Array(0.33333333, 0.33333333, 0.33333333, 0.05961587, 0.47014206, 0.47014206, 0.47014206, 0.05961587, 0.47014206, 0.47014206, 0.47014206, 0.05961587, 0.79742699, 0.10128651, 0.10128651, 0.10128651, 0.79742699, 0.10128651, 0.10128651, 0.10128651, 0.79742699)
        
        Dim quadPt(2), coord
        Dim area: area = TriangleArea(arrVerts)
        Dim arrQuadratureForce, arrQuadratureMoment, arrQuadratureRadius, arrSum
        Dim quadID: quadID = 0
        arrSum = Array(0, 0, 0)
        While quadID <= UBound(quadWts)
                For coord= 0 To 2
                        quadPt(coord) = (arrVerts(0)(coord) * quadCoords(quadID * 3)) + (arrVerts(1)(coord) * quadCoords(quadID * 3 + 1)) + (arrVerts(2)(coord) * quadCoords(quadID * 3 + 2))
                Next
                'retrieve the component stresses at the quadrature point
                arrQuadratureForce = Array(objSnSPlugin.QuerySolutionValue(strObject, "SIGXX", quadPt), objSnSPlugin.QuerySolutionValue(strObject, "SIGYY", quadPt), objSnSPlugin.QuerySolutionValue(strObject, "SIGZZ", quadPt))
                'multiply component stress by area to get force at quadrature point
                arrQuadratureForce = Rhino.VectorScale(arrQuadratureForce, quadWts(quadID) * area)
                'compute the radius to the quadrature point
                arrQuadratureRadius = Rhino.VectorSubtract(quadPt, arrCentroid)
                'compute RxF to get the contribution to the moment
                arrQuadratureMoment = Rhino.VectorCrossProduct(arrQuadratureRadius, arrQuadratureForce)
                'accumulate the sum for the numerical integration, skip numerical outliers
                If Rhino.VectorLength(arrQuadratureForce) < 10000 Then
                        arrSum = Rhino.VectorAdd(arrSum, arrQuadratureMoment)
                End If
                
                quadID = quadID + 1
        Wend
        IntegrateMoments = arrSum
End Function

'TriangleArea: Utility function to compute the area of the triangle defined 
'by the array of vertices in arrVerts
Function TriangleArea(arrVerts)
        Dim vec0: vec0 = Rhino.VectorSubtract(arrVerts(1), arrVerts(0))
        Dim vec1: vec1 = Rhino.VectorSubtract(arrVerts(2), arrVerts(0))
        Dim vecCross: vecCross = Rhino.VectorCrossProduct(vec0, vec1)
        TriangleArea = 0.5 * Rhino.VectorLength(vecCross)
End Function

'StressForce: Utility function to integrate a stress components over a surface
'to determine the moments acting about the centroid of the surface
'strObject -the solid with SnS solution data
'strSurfaceObject -the surface over which the integration should occur (typically a cross-section of strObject)
Function StressMoment(strSolid, strSection)
        Dim meshObject,arrFaces,arrVerts(2),fid,arrCentroid
        Dim sumStressMoment
        sumStressMoment = Array(0, 0, 0)
        
        'Extract the analysis mesh for the surface objects
        If Rhino.IsObject(strSection) Then
                meshObject = Rhino.ExtractAnalysisMesh(strSection)
        Else
                Rhino.Print "Invalid surface object"
                StressMoment = Array(0, 0, 0)
        End If

        'Integrate the desired stress component over the facets of the extracted mesh
        If Not IsNull(meshObject) Then
                arrCentroid = Rhino.SurfaceAreaCentroid(strSection)
                Rhino.MeshQuadsToTriangles meshObject
                arrFaces = Rhino.MeshFaces(meshObject, False)
                If IsArray(arrFaces) Then
                        fid = 0
                        While fid <= UBound(arrFaces)
                                arrVerts(0) = arrFaces(fid)
                                arrVerts(1) = arrFaces(fid + 1)
                                arrVerts(2) = arrFaces(fid + 2)
                                fid = fid + 3
                                sumStressMoment = Rhino.VectorAdd(sumStressMoment, IntegrateMoments(strSolid, arrVerts, arrCentroid(0)))
                        Wend
                End If
                Rhino.DeleteObject meshObject
        Else
                Rhino.Print "Unable to extract analysis mesh"
        End If
        StressMoment = sumStressMoment
End Function

Call Main()
Sub Main()
        Dim strObject, strPlane, arrCurves, arrSection, strSection, arrMoments
        On Error Resume Next
        
        Set objSnSPlugIn = Rhino.GetPluginObject("SnSScript")
        objSnSPlugIn.SetOutputOnOff False
        
        If Err Then 
                MsgBox Err.Description
                Exit Sub
        End If
        
        strObject = Rhino.GetObject("Select a solid")
        strPlane = Rhino.GetObject("Select a plane", 8)
        
        'compute the curves of the intersection between the solid and the plane
        arrCurves = Rhino.IntersectBreps(strObject, strPlane)

        If Not IsNull(arrCurves) Then
                'construct a planar surface
                arrSection = Rhino.AddPlanarSrf(arrCurves)
                'WARNING!  The output from AddPlanarSrf can be a collection of surfaces.
                'For simplicity only the zeroth one is used.  Otherwise their collective
                'centroid must be used and their results accumulated.  This code DOESN'T
                'do that.
                arrMoments = StressMoment(strObject, arrSection(0))
                Rhino.Print arrMoments(0) & "," & arrMoments(1) & "," & arrMoments(2)
                'clean up by deleting curves and the cross-sections
                Dim i
                For i = 0 To ubound(arrCurves)
                        Rhino.DeleteObject arrCurves(i)
                Next
                For i= 0 To ubound(arrSection)
                        Rhino.DeleteObject arrSection(i)
                Next
        End If
        
End Sub

Line #	Description
1	Only scripts with Option Explicit specified can be debugged in the RhinoScript editor (Monkey).
6	Global variable for the SnSScript plug-in object.
13-42	Define a function to numerically integrate a solution component over a triangle.
14-15	Define the quadrature weights and points for integrating over a triangle.
18	Compute the area of the current triangle.
22-40	Iterate over quadrature weights and points, evaluating the integral of the bending moment due to component stresses.
46-51	Define a function to compute the area of a triangle using the cross-product.
57-90	Define a function to integrate moments due to component stresses over a surface in Rhino.
64	Extract a mesh approximation of the surface.
73	Convert the mesh quads to triangles for compatibility with integration algorithm.
77-83	Integrate over each mesh face, accumulating their contributions in sumStressMoment.
85	Delete the mesh object for the surface.
93-130	The main body of the script.
97	Load the SnSScript plug-in object.
105	Prompt the user to select a solid from the Rhino document. The solid is identified by its GUID stored in strObject.
106	Prompt the user to select a section plane from the Rhino document.
109	Intersect the plane with the solid to get the intersection curves.
113	Convert the section curves into a trimmed plane.
118-119	Compute the bending moments from the component stresses and print them out.
121-127	Clean up the temporary objects.