**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