Tetrahedra discretization of a beam

In this post I will introduce the advances of the last weeks regarding the discretization into tetrahedra of a beam with Grasshopper's VB component. Rotations about an arbitrary axis An important part of the research has been devoted to the topic of 3d rotation, as for locating the sub-nodes of the beam appropiately we need to transform the beam's local coordinates with respect of its axis into the actual world's coordinates for Rhino to interpret.
The most popular way to make it are rotation matrices (here there is a thorough explanation on how), but also quaternions have been contemplated and still not fully discarded. Apparently rotation matrices present some numerical drifting problems when utilized for many time steps, as well as other problems that quaternions solve quite easily (mentioned here). Tetrahedra subdivision The beam is defined by two points arbitrarily located in 3d space. These points define a director vector and also a distance. The remaining geometric properties (height h and width b) in the attached code have been assimilated to a percentage of the lenght. Also the number of segments s has been hardcoded to be proportional to the width but it could be anything else. The attached code uses two functions: one for coordinate transformation by means of transformation matrices and another for matrix multiplication. The rest of the subroutine simply locates points along the beam according to the width b, the height h and the segment separation s.
The image above shows how the beam is discretized into nodes located at the edges. Each segment gives place to five embedded tetrahedron (right). Local axes of the beam will be later utilized for the matter integration. Code
Private Sub RunScript(ByVal ni As Point3d, ByVal nj As Point3d, ByVal angle As Double, ByRef geom As Object)
    'This code was written by Rabindranath Andujar, Jordi Casabo, Jaume Roset, Vojko Kilar and Simon Petrovcic with the purposes 'of the development of Rabindranath Andujar's Phd Thesis "Stochastic Simulation and Lagrangian dynamics applied to Structural 'Design" 'Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, 'Version 1.3 Or any later version published by the Free Software Foundation; With no Invariant Sections, no Front-Cover Texts, 'And no Back-Cover Texts. 'This code is intended to facilitate a tetrahedric discretization of a beam provided two points (initial and final). It should 'serve to be inserted into a wider matter integration code (FEM/FVM/MSS, etc) for which nodal coordinates, topological table 'and material properties are eventually given. Dim L As Double 'lenght of the beam Dim s As Double 'space between divisions Dim b As Double 'width of the beam Dim h As Double 'height of the beam Dim u,v,w As Double 'components of the director vector of the beam Dim lines As New list(Of rhino.geometry.line) 'geometric representation of the beam Dim Nodes As New List(Of rhino.geometry.point3d) 'geometric nodes of the beam, where matter will be integrated Dim n As Integer 'number of nodes Dim i As Integer 'auxiliar index Dim line As Rhino.Geometry.Line 'auxiliar instantiation variable Dim Node As rhino.geometry.point3d 'Auxiliar coordinates for the shape of the beam 'The components of the director vector of the beam are the difference between the initial and the final points u = nj.X - ni.X v = nj.y - ni.y w = nj.z - ni.z ' The lenght of the beam is the norm of the vector L = math.Sqrt(u * u + v * v + w * w) 'We have chosen arbitrarily here a height h=10% of the lenght, width b=60% of the height and spacing a portion of the width h = L * 0.1 b = h * 0.6 n = int(L / b) s = L / n 'The nodes located in the vertices of the section are positioned along the beam by means of matrix 'transformations with respect to the director vector and the initial point of the beam. 'The function TransformCoordinates serves such purpose. 'u,v,w coordinates enter divided by the lenght of the beam in order to operate with a normalized vector. 'The angle parameter represents the rotation of the section around the longitudinal axis (local x) ' local z axis '    | '    |   angle of rotation along x local axis '    |  / '  *---* ' h| |/| ' -|-+-|---local y axis '  | | | '  *---* '    b For i = 0 To n 'Top left corner of the beam section Node = TransformCoordinates(ni.X, ni.Y, ni.Z, -b / 2, h / 2, S * i, u / L, v / L, w / L, angle) 'Then we collect their values into Rhino Point3Ds Nodes.Add(Node) 'Top right corner of the beam section Node = TransformCoordinates(ni.X, ni.Y, ni.Z, b / 2, h / 2, S * i, u / L, v / L, w / L, angle) 'Then we collect their values into Rhino Point3Ds Nodes.Add(Node) 'Lower right corner of the beam section Node = TransformCoordinates(ni.X, ni.Y, ni.Z, b / 2, -h / 2, S * i, u / L, v / L, w / L, angle) 'Then we collect their values into Rhino Point3Ds Nodes.Add(Node) 'Lower left corner of the beam section Node = TransformCoordinates(ni.X, ni.Y, ni.Z, -b / 2, -h / 2, S * i, u / L, v / L, w / L, angle) 'Then we collect their values into Rhino Point3Ds Nodes.Add(Node) Next i 'This loop repeats the cross section of the beam (a b x h rectangle) using the vertices previously collected 'along a beam defined by an initial and a final points. 'Also the longitudinal edges of the prism are iterated for every segment of beam. For i = 0 To nodes.Count - 8 Step 4 'Four sides of the section Line = New Rhino.Geometry.Line(nodes(i), nodes(i + 1)) Lines.Add(Line) Line = New Rhino.Geometry.Line(nodes(i + 1), nodes(i + 2)) Lines.Add(Line) Line = New Rhino.Geometry.Line(nodes(i + 2), nodes(i + 3)) Lines.Add(Line) Line = New Rhino.Geometry.Line(nodes(i + 3), nodes(i)) Lines.Add(Line) 'Longitudinal segments of the edges Line = New Rhino.Geometry.Line(nodes(i), nodes(i + 4)) Lines.Add(Line) Line = New Rhino.Geometry.Line(nodes(i + 1), nodes(i + 5)) Lines.Add(Line) Line = New Rhino.Geometry.Line(nodes(i + 2), nodes(i + 6)) Lines.Add(Line) Line = New Rhino.Geometry.Line(nodes(i + 3), nodes(i + 7)) Lines.Add(Line) Next i 'Here we iterate in order to make the outer longitudinal diagonals and the inner cross sectional ones 'On each segment we have to change the direction, so they are alternated. It was carefully done, as the 'final purpose is to create the edges of tetrahedra that do not overlap. 'To such end, diagonals change every two segments in outer faces as in inner cross sectional connections. For i = 1 To nodes.Count - 4 Step 8 'Top faces Line = New Rhino.Geometry.Line(nodes(i), nodes(i + 5)) Lines.Add(Line) Line = New Rhino.Geometry.Line(nodes(i + 5), nodes(i + 8)) Lines.Add(Line) 'Right side faces Line = New Rhino.Geometry.Line(nodes(i + 2), nodes(i + 5)) Lines.Add(Line) Line = New Rhino.Geometry.Line(nodes(i + 5), nodes(i + 10)) Lines.Add(Line) 'Bottom faces Line = New Rhino.Geometry.Line(nodes(i + 2), nodes(i + 3)) Lines.Add(Line) Line = New Rhino.Geometry.Line(nodes(i + 3), nodes(i + 10)) Lines.Add(Line) 'Left side faces Line = New Rhino.Geometry.Line(nodes(i), nodes(i + 3)) Lines.Add(Line) Line = New Rhino.Geometry.Line(nodes(i + 3), nodes(i + 8)) Lines.Add(Line) 'Sectional diagonals Line = New Rhino.Geometry.Line(nodes(i), nodes(i + 2)) 'Lines.Add(Line) Line = New Rhino.Geometry.Line(nodes(i + 3), nodes(i + 5)) 'Lines.Add(Line) Next i 'And the final section to close the prism: perimetral and diagonal lines Line = New Rhino.Geometry.Line(nodes(nodes.count - 4), nodes(nodes.count - 3)) Lines.Add(Line) Line = New Rhino.Geometry.Line(nodes(nodes.count - 3), nodes(nodes.count - 2)) Lines.Add(Line) Line = New Rhino.Geometry.Line(nodes(nodes.count - 2), nodes(nodes.count - 1)) Lines.Add(Line) Line = New Rhino.Geometry.Line(nodes(nodes.count - 1), nodes(nodes.count - 4)) Lines.Add(Line) 'Diagonal goes from the last node to the node in the opposite side of the section Line = New Rhino.Geometry.Line(nodes(nodes.count - 1), nodes(nodes.count - 3)) Lines.Add(Line) 'Finally everything is delivered to Rhino for visualization Geom = Lines End Sub
  Function TransformCoordinates(a As Double, b As Double, c As Double, x As Double, y As Double, z As Double, u As Double, v As Double, w As Double, Alpha As Double) As rhino.Geometry.Point3d 'This function returns a Rhino.Geometry.Point3d rotated around an arbitrary axis. 'The initial point of the beam is defined by a,b and c coordinates 'Point's coordinates to be transformed are x,y,z 'The axis is defined by the u,v,w coordinates, which should be normalized 'Angle is provided in the Alpha parameter Dim CAlpha As Double 'Auxiliar variable for storing the cosine of the rotation around the axial axis Dim SAlpha As Double 'Auxiliar variable for storing the sine of the rotation around the axial axis Dim d As Double 'Auxiliar variable for storing the sine of the rotation around the y axis Dim Vector(3) As Double 'Auxiliar Vector to store and make the neccessary operations Dim Point As rhino.Geometry.Point3d 'The returned instance of a modified Point3d 'Initialize variables Vector(0) = x Vector(1) = y Vector(2) = z Vector(3) = 1 d = math.Sqrt(v ^ 2 + w ^ 2) CAlpha = math.Cos(Alpha) SAlpha = math.Sin(Alpha) 'The translation matrix moves the point to the a,b,c coordinates Dim Transl(,) As Double = { _ {1, 0, 0, a}, _ {0, 1, 0, b}, _ {0, 0, 1, c}, _ {0, 0, 0, 1}} 'Rotation matrix along the world's x axis Dim Rotx(,) As Double = { _ {1, 0, 0, 0}, _ {0, w/d, v/d, 0}, _ {0, -v/d, w/d, 0}, _ {0, 0, 0, 1}} 'Rotation matrix along the beam's y axis Dim Roty(,) As Double = { _ {d, 0, u, 0}, _ {0,1, 0, 0}, _ {-u, 0, d, 0}, _ {0, 0, 0, 1}} 'Rotation matrix along the world's z axis Dim Rotz(,) As Double = { _ {CAlpha, -SAlpha, 0, 0}, _ {SAlpha, CAlpha, 0, 0}, _ {0, 0, 1, 0}, _ {0, 0, 0, 1}} 'Each transformation is made sequentially as a matrix-vector product Vector = MultiplyMatrixVector(Rotz, Vector) Vector = MultiplyMatrixVector(Roty, Vector) Vector = MultiplyMatrixVector(Rotx, Vector) Vector = MultiplyMatrixVector(Transl, Vector) 'And eventually data is returned into the auxiliar point3d Point.X = Vector(0) Point.y = Vector(1) Point.z = Vector(2) Return point End Function
  Function MultiplyMatrixVector(Matrix(,) As Double, Vector() As Double) As Double()     'This function returns the product of a matrix times a vector.     'Matrix column number must be the same as vector's components number:     '     '                | |             | |     '  |       |     | |             | |     'M=|  mxn  | , V=|n|   ===>  M·V=|n|     '  |       |     | |             | |     '                | |             | |     '     Dim i,j As Integer 'Auxiliar indexes     Dim n As Integer 'Matrix's column number     Dim VAux() As Double 'Auxiliar temporary vector     'Vector's size is reset to the number of matrix's columns     n = Ubound(Matrix, 1)     ReDim VAux(n)     'In this loop we proceed to the summation of the products of the matrix terms times the     'corresponding vector terms.     'First we iterate by the columns     For i = 0 To n       'And then by each term of the vector       For j = 0 To n         'Auxiliar vector accumulates the products in the corresponding index         Vaux(i) = Vaux(i) + Matrix(i, j) * Vector(j)       Next j     Next     'And then we return the resulting vector     Return Vaux   End Function