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