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

Yet one more article...

And another article published!
This once in the journal Iranian journal WASJ (World Applied Science Journal).
This is the link: www.idosi.org/wasj/wasj14(8)11/21.pdf.

Abstract: Structural dynamics is a rather complex field of research that concerns to a broad range of disciplines, from structural engineering to graphics animation, robotics or aeronautics. A primary consequence of this is an overwhelming amount of literature on the topic, apparently disconnected, as each author focuses on his / her particular field. To complicate things further, the daunting list of numerical methods severely blurs the scope of the researcher, making it very difficult to understand what their purpose is in each case and even if these are applicable to the analysis of structural behavior. This paper presents a reference framework where researchers and developers from diverse disciplines can assess the main methods currently used in structural dynamics simulation. A direct correlation is made between methods to solve  Ordinary, Partial and Algebraic Differential Equations and their physical counterparts Time, Matter and Constraints. It is also discussed their application in different industries.



Thanks to my mentors, Jaume and Vojko
 

My first SCI Indexed Article!

Here it is! A whole achievement!

After some while working hard in collaboration with the Faculty of Architecture of  the University of Ljubljana the result of our research appears published in a scientific Journal: http://www.ttem-bih.org/ttem_3_2011.html.
As a computational dynamics rookie myself I find it very hard to understand the whole complexity of the topic. For that reason, we decided to make an introductory overview of what is needed to start scouting through such a steepy and accidented path. I hope it serves to somebody else...

However, from all the article my favourite part is at the beginning: 

“Structural engineering is the art of molding materials we don’t wholly understand, into shapes we can’t fully analyze, so as to withstand forces we can’t really assess, in such a way that the community at large has no reason to suspect the extent of our ignorance.”

Which I found as a quote from James E. Amrhein of the Masonry Institute of America...

Well, happy autumn!

First FEM + Lagrange multipliers program in Grasshopper

Good, good...summer is about to finish and here I am, back in Ljubljana!

The last post was created two months ago and a shift in the researching techinques was announced: from C++ under Code:Blocks, Visual Studio or any other IDE (never found the one that fitted my needs), to Visual Basic under Grasshopper...

What a change! It really has boosted productivity...debugging is now as easy as dragging a panel and linking it to a component, results of the changes in the program can be viewed on real time and also with a bit of tweaking its consequences in a CAD model can be represented...wow..absolutely good...

In concrete means, during the summer I have been able to produce the whole Grasshopper definition shown below:

The initial part simply parses data from three files that I chose to be in .txt format, but anyone else could have chosen any other way (i.e. Excel or OpenOffice): Constraints.txt, Nodes.txt and ExternalForces.txt.
The heavy programming components are FEM, LM and LUPD, and I provide their contents here:

FEM Component:
Private Sub RunScript(ByVal Constraints As List(Of String), ByVal Nodes As List(Of String), ByVal Forces As List(Of String), ByRef Kg As Object, ByRef u As Object, ByRef F As Object, ByRef Geom As Object) 'your code goes here… Dim Lenght As Double Dim Lxy As Double Dim ni As Int16 Dim nj As int16 Dim EModulus As Double Dim GModulus As Double Dim Area As Double Dim Theta As Double Dim InertiaY As Double Dim InertiaZ As Double Dim JModulus As Double Dim PosNi(2) As Double Dim PosNj(2) As Double Dim i As Int16 Dim j As Int16 Dim k As Int16 Dim m As int16 Dim n As Int16 Dim ConstraintType As int32 Dim lines As New list(Of rhino.geometry.line) Dim KgAux(,) As Double ReDim KgAux(6 * (nodes.count - 1), 6 * (nodes.count - 1)) For i = 1 To Constraints.Count - 1 ConstraintType = CInt(constraints(i).split(",")(0)) If ConstraintType = 0 Then ni = CInt(constraints(i).split(",")(1)) nj = CInt(constraints(i).split(",")(2)) Emodulus = CDbl(constraints(i).split(",")(3)) GModulus = CDbl(constraints(i).split(",")(4)) Area = CDbl(constraints(i).split(",")(5)) InertiaY = CDbl(constraints(i).split(",")(6)) InertiaZ = CDbl(constraints(i).split(",")(7)) JModulus = CDbl(constraints(i).split(",")(8)) Theta = CDbl(constraints(i).split(",")(9)) posni(0) = CDbl(nodes(ni).split(",")(0)) posni(1) = CDbl(nodes(ni).split(",")(1)) posni(2) = CDbl(nodes(ni).split(",")(2)) posnj(0) = CDbl(nodes(nj).split(",")(0)) posnj(1) = CDbl(nodes(nj).split(",")(1)) posnj(2) = CDbl(nodes(nj).split(",")(2)) Dim Line As New Rhino.Geometry.Line(posni(0), posni(1), posni(2), posnj(0), posnj(1), posnj(2)) Lines.Add(Line) Lenght = math.Sqrt((posnj(0) - posni(0)) ^ 2 + (posnj(1) - posni(1)) ^ 2 + (posnj(2) - posni(2)) ^ 2) Lxy = math.Sqrt((posnj(0) - posni(0)) ^ 2 + (posnj(1) - posni(1)) ^ 2) 'Lets compose the element's matrix first Dim Ke(11,11) As Double For j = 0 To 5 For k = 0 To 5 Ke(j, k) = kenini(Lenght, Emodulus, GModulus, Area, InertiaY, InertiaZ, JModulus)(j, k) Ke(j, 6 + k) = keninj(Lenght, Emodulus, GModulus, Area, InertiaY, InertiaZ, JModulus)(j, k) Ke(6 + j, k) = kenjni(Lenght, Emodulus, GModulus, Area, InertiaY, InertiaZ, JModulus)(j, k) Ke(6 + j, 6 + k) = kenjnj(Lenght, Emodulus, GModulus, Area, InertiaY, InertiaZ, JModulus)(j, k) Next k Next j 'Then create the transformation matrix components Dim s1,s2,s3 As Double Dim c1,c2,c3 As Double Dim Tinv(11,11) As Double Dim T(11,11) As Double If Lxy > 0.0001 Then s1 = (posnj(1) - posni(1)) / Lxy c1 = (posnj(0) - posni(0)) / Lxy Else s1 = 0.0 c1 = 1.0 End If s2 = (posnj(2) - posni(2)) / Lenght c2 = Lxy / Lenght s3 = math.sin(theta) c3 = math.cos(theta) T = Transformatrix(c1 * c2, s1 * c2, s2, -c1 * s2 * s3 - s1 * c3, -s1 * s2 * s3 + c1 * c3, s3 * c2, -c1 * s2 * c3 - s1 * s3, -s1 * s2 * c3 - c1 * s3, c3 * c2) Tinv = TransformatrixT(c1 * c2, s1 * c2, s2, -c1 * s2 * s3 - s1 * c3, -s1 * s2 * s3 + c1 * c3, s3 * c2, -c1 * s2 * c3 - s1 * s3, -s1 * s2 * c3 - c1 * s3, c3 * c2) 'And proceed to transform the element's matrix into global coordinates Dim KeAux(11,11) As Double KeAux = multiply2DMatrices(Ke, Tinv) Ke = multiply2DMatrices(T, KeAux) For j = 0 To 11 print(ke(j, 0) & "," & ke(j, 1) & "," & ke(j, 2) & "," & ke(j, 3) & "," & ke(j, 4) & "," & ke(j, 5) & "," & ke(j, 6) & "," & ke(j, 7) & "," & ke(j, 8) & "," & ke(j, 9) & "," & ke(j, 10) & "," & ke(j, 11)) Next print("****************") 'Once translated into global coordinates, we can assemble the element in the temporary global matrix For j = 0 To 5 For k = 0 To 5 KgAux(6 * (ni - 1) + j, 6 * (ni - 1) + k) = KgAux(6 * (ni - 1) + j, 6 * (ni - 1) + k) + Ke(j, k) KgAux(6 * (ni - 1) + j, 6 * (nj - 1) + k) = KgAux(6 * (ni - 1) + j, 6 * (nj - 1) + k) + Ke(j, 6 + k) KgAux(6 * (nj - 1) + j, 6 * (ni - 1) + k) = KgAux(6 * (nj - 1) + j, 6 * (ni - 1) + k) + Ke(6 + j, k) KgAux(6 * (nj - 1) + j, 6 * (nj - 1) + k) = KgAux(6 * (nj - 1) + j, 6 * (nj - 1) + k) + Ke(6 + j, 6 + k) Next k Next j End If Next i 'And parse it into a tree structure for management into Grasshopper Dim Tree As New DataTree(Of Double) For m = 0 To ubound(KgAux) - 1 'Iterate over the second dimension of the array For n = 0 To ubound(KgAux) - 1 'Insert the current value into the DataTree. 'The first dimension ends up as the Path, 'the second dimension ends up as the index. tree.Insert(KgAux(m, n), New GH_Path(m), n) Next n Next m Kg = tree F = ReadForces(Forces) u = InitializeDisplacements(Nodes) Geom = Lines End Sub '<Custom additional code> Function KeNiNi(L As Double, E As Double, G As Double, A As Double, Iy As Double, Iz As Double, J As Double) As Double(,) Dim Kaux(,) As Double = { _ {A * E / L, 0, 0, 0, 0, 0}, _ {0, 12 * E * Iz / (L * L * L), 0, 0, 0, 6*E*Iz/(L * L)}, _ {0, 0, 12 * E * Iy / (L * L * L), 0, -6*E*Iy/(L * L), 0}, _ {0, 0, 0, G * J / L, 0, 0}, _ {0, 0, -6*E*Iy/(L * L), 0, 4*E*Iy/L, 0}, _ {0, 6*E*Iz/(L * L), 0, 0, 0, 4*E*Iz/L} _ } Return Kaux End Function Function KeNiNj(L As Double, E As Double, G As Double, A As Double, Iy As Double, Iz As Double, J As Double) As Double(,) Dim Kaux(,) As Double = { _ {-A*E/L, 0, 0, 0, 0, 0}, _ {0, -12*E*Iy/(L*L*L), 0, 0, 0, 6*E*Iz/(L*L)}, _ {0, 0, -12*E*Iz/(L*L*L), 0, -6*E*Iy/(L*L),0}, _ {0, 0, 0,-G*J/L, 0, 0}, _ {0, 0, 6*E*Iy/(L*L), 0, 2*E*Iy/L, 0}, _ {0, -6*E*Iz/(L*L), 0, 0, 0, 2*E*Iz/L} _ } Return Kaux End Function Function KeNjNi(L As Double, E As Double, G As Double, A As Double, Iy As Double, Iz As Double, J As Double) As Double(,) Dim Kaux(,) As Double = { _ {-A*E/L, 0, 0, 0, 0, 0}, _ {0, -12 * E * Iz / (L * L * L), 0, 0, 0, -6*E*Iz/(L * L)}, _ {0, 0, -12 * E * Iy / (L * L * L), 0, 6*E*Iy/(L * L), 0}, _ {0, 0, 0, -G * J / L, 0, 0}, _ {0, 0, -6*E*Iy/(L * L), 0, 2*E*Iy/L, 0}, _ {0, 6*E*Iz/(L * L), 0, 0, 0, 2*E*Iz/L} _ } Return Kaux End Function Function KeNjNj(L As Double, E As Double, G As Double, A As Double, Iy As Double, Iz As Double, J As Double) As Double(,) Dim Kaux(,) As Double = { _ {A*E/L, 0, 0, 0, 0, 0}, _ {0, 12 * E * Iz / (L * L * L), 0, 0, 0, -6*E*Iz/(L * L)}, _ {0, 0, 12 * E * Iy / (L * L * L), 0, 6*E*Iy/(L * L), 0}, _ {0, 0, 0, G * J / L, 0, 0}, _ {0, 0, 6*E*Iy/(L * L), 0, 4*E*Iy/L, 0}, _ {0, -6*E*Iz/(L * L), 0, 0, 0, 4*E*Iz/L} _ } Return Kaux End Function Function Transformatrix(lx As Double, mx As Double, nx As Double, ly As Double, my As Double, ny As Double, lz As Double, mz As Double, nz As Double) As Double(,) Dim Kaux(,) As Double = { _ {lx, mx, nx, 0, 0, 0, 0, 0, 0, 0, 0, 0}, _ {ly, my, ny, 0, 0, 0, 0, 0, 0, 0, 0, 0}, _ {lz, mz, nz, 0, 0, 0, 0, 0, 0, 0, 0, 0}, _ {0, 0, 0, lx, mx, nx, 0, 0, 0, 0, 0, 0}, _ {0, 0, 0, ly, my, ny, 0, 0, 0, 0, 0, 0}, _ {0, 0, 0, lz, mz, nz, 0, 0, 0, 0, 0, 0}, _ {0, 0, 0, 0, 0, 0, lx, mx, nx, 0, 0, 0}, _ {0, 0, 0, 0, 0, 0, ly, my, ny, 0, 0, 0}, _ {0, 0, 0, 0, 0, 0, lz, mz, nz, 0, 0, 0}, _ {0, 0, 0, 0, 0, 0, 0, 0, 0, lx, mx, nx}, _ {0, 0, 0, 0, 0, 0, 0, 0, 0, ly, my, ny}, _ {0, 0, 0, 0, 0, 0, 0, 0, 0, lz, mz, nz} _ } Return Kaux End Function Function TransformatrixT(lx As Double, mx As Double, nx As Double, ly As Double, my As Double, ny As Double, lz As Double, mz As Double, nz As Double) As Double(,) Dim Kaux(,) As Double = { _ {lx, ly, lz, 0, 0, 0, 0, 0, 0, 0, 0, 0}, _ {mx, my, mz, 0, 0, 0, 0, 0, 0, 0, 0, 0}, _ {nx, ny, nz, 0, 0, 0, 0, 0, 0, 0, 0, 0}, _ {0, 0, 0, lx, ly, lz, 0, 0, 0, 0, 0, 0}, _ {0, 0, 0, mx, my, mz, 0, 0, 0, 0, 0, 0}, _ {0, 0, 0, nx, ny, nz, 0, 0, 0, 0, 0, 0}, _ {0, 0, 0, 0, 0, 0, lx, ly, lz, 0, 0, 0}, _ {0, 0, 0, 0, 0, 0, mx, my, mz, 0, 0, 0}, _ {0, 0, 0, 0, 0, 0, nx, ny, nz, 0, 0, 0}, _ {0, 0, 0, 0, 0, 0, 0, 0, 0, lx, ly, lz}, _ {0, 0, 0, 0, 0, 0, 0, 0, 0, mx, my, mz}, _ {0, 0, 0, 0, 0, 0, 0, 0, 0, nx, ny, nz} _ } Return Kaux End Function Function Multiply2DMatrices(A(,) As Double, B(,) As Double) As Double(,) Dim i As int32 Dim j As int32 Dim k As int32 Dim MAux(,) As Double ReDim MAux(ubound(A, 1), ubound(B, 2)) For i = 0 To ubound(A, 1) For j = 0 To ubound(B, 1) For k = 0 To ubound(B, 2) MAux(i, j) = MAux(i, j) + A(i, k) * B(k, j) Next k Next Next Return MAux End Function Function ReadForces(F As list(Of String)) As datatree(Of Double) Dim i As Integer Dim FAux As New datatree(Of Double) For i = 1 To F.count - 1 FAux.insert(CDbl(F(i).split(",")(1)), New GH_Path(i), 0) FAux.insert(CDbl(F(i).split(",")(2)), New GH_Path(i), 1) FAux.insert(CDbl(F(i).split(",")(3)), New GH_Path(i), 2) FAux.insert(CDbl(F(i).split(",")(4)), New GH_Path(i), 3) FAux.insert(CDbl(F(i).split(",")(5)), New GH_Path(i), 4) FAux.insert(CDbl(F(i).split(",")(6)), New GH_Path(i), 5) Next Return FAux End Function Function InitializeDisplacements(N As list(Of String)) As datatree(Of Double) Dim i As Integer Dim NAux As New datatree(Of Double) For i = 1 To N.count - 1 NAux.insert(0, New GH_Path(i), 0) NAux.insert(0, New GH_Path(i), 1) NAux.insert(0, New GH_Path(i), 2) NAux.insert(0, New GH_Path(i), 3) NAux.insert(0, New GH_Path(i), 4) NAux.insert(0, New GH_Path(i), 5) Next Return NAux End Function

LM Component:
Private Sub RunScript(ByVal Constraints As List(Of String), ByVal Kg As DataTree(Of Double), ByVal u As DataTree(Of Double), ByVal F As DataTree(Of Double), ByRef KgExt As Object, ByRef uExt As Object, ByRef FExt As Object) 'your code goes here… Dim ConstraintType As Int16 Dim RotationX As Double Dim RotationY As Double Dim RotationZ As Double Dim DisplacementX As Double Dim DisplacementY As Double Dim DisplacementZ As Double Dim i As int32 Dim j As int32 Dim k As int32 Dim ExtraRow As int32 Dim ExtraCol As int32 Dim ni As int32 Dim KgAux As New DataTree(Of Double)(Kg) For i = 1 To Constraints.Count - 1 ConstraintType = CInt(constraints(i).split(",")(0)) ni = CInt(constraints(i).split(",")(1)) If ConstraintType = 1 Then RotationX = CDbl(constraints(i).split(",")(3)) RotationY = CDbl(constraints(i).split(",")(4)) RotationZ = CDbl(constraints(i).split(",")(5)) DisplacementX = CDbl(constraints(i).split(",")(6)) DisplacementY = CDbl(constraints(i).split(",")(7)) DisplacementZ = CDbl(constraints(i).split(",")(8)) InsertConstraint(DisplacementX, KgAux, 6 * (ni - 1)) InsertConstraint(DisplacementY, KgAux, 6 * (ni - 1) + 1) InsertConstraint(DisplacementZ, KgAux, 6 * (ni - 1) + 2) InsertConstraint(RotationX, KgAux, 6 * (ni - 1) + 3) InsertConstraint(RotationY, KgAux, 6 * (ni - 1) + 4) InsertConstraint(RotationZ, KgAux, 6 * (ni - 1) + 5) End If Next i ExtraCol = KgAux.BranchCount ExtraRow = KgAux.BranchCount For i = Kg.BranchCount To ExtraRow - 1 For j = Kg.BranchCount To ExtraCol - 1 KgAux.Insert(0, New GH_Path(i), j) Next j Next i j = F.BranchCount + 1 k = 0 FExt = New DataTree(Of Double)(F) For i = Kg.BranchCount To ExtraRow - 1 FExt.insert(0, New GH_Path(j), k) k = k + 1 Next uExt = New DataTree(Of Double)(u) KgExt = New DataTree(Of Double)(KgAux) 'Dim Txt As New System.IO.StreamWriter("C:\Users\John\Desktop\Matrix1.csv") 'For i = 0 To KgAux.BranchCount - 1 'For j = 0 To KgAux.BranchCount - 1 ' print(i & "," & j) ' Txt.Write(KgAux.Branch(i).item(j) & ",") ' Next ' txt.WriteLine() 'Next 'txt.Flush() End Sub '<Custom additional code> Sub InsertConstraint(Value As Double, ByRef K As datatree(Of Double), n As int32) Dim i As int32 Dim ExtraCol As int32 Dim ExtraRow As int32 ExtraCol = K.BranchCount ExtraRow = K.BranchCount For i = 0 To n - 1 K.Insert(0, New GH_Path(ExtraRow), i) K.Insert(0, New GH_Path(i), ExtraCol) Next i K.Insert(Value, New GH_Path(ExtraRow), n) K.Insert(Value, New GH_Path(n), ExtraCol) For i = n + 1 To K.BranchCount - 1 K.Insert(0, New GH_Path(ExtraRow), i) K.Insert(0, New GH_Path(i), ExtraCol) Next i End Sub

LUPD Component:
Private Sub RunScript(ByVal A As DataTree(Of Double), ByVal X As DataTree(Of Double), ByVal B As DataTree(Of Double), ByRef Xsol As Object) 'your code goes here… Dim n As Int32 Dim er As Int32 Dim i As Int32 Dim j As Int32 Dim m As int32 Dim BAux() As Double Dim XAux() As Double Dim tol As Double Dim AAux(,) As Double n = A.BranchCount - 1 ReDim Aaux(n, n) ReDim Baux(n) ReDim Xaux(n) For i = 0 To n For j = 0 To n AAux(i, j) = A.Branch(i).item(j) Next j Next i For i = 0 To B.BranchCount - 1 For j = 0 To b.Branch(i).Count - 1 BAux(m) = B.Branch(i).item(j) m = m + 1 Next j Next i tol = 1E-17 LUD(AAux, BAux, n, XAux, tol, er) Xsol = XAux End Sub '<Custom additional code> Sub LUD(a(,) As Double, b() As Double, n As Integer, x() As Double, tol As Double, er As Int32) Dim o() As Double Dim s() As Double ReDim o(n) ReDim s(n) Call Decompose(a, n, tol, o, s, er) If er = 0 Then Call Substitute(a, o, n, b, x) Else MsgBox("ill-conditioned system") Exit Sub End If End Sub Sub Decompose(a(,) As Double, n As int32, tol As Double, o() As Double, s() As Double, er As int32) Dim i As Integer Dim j As Integer Dim k As Integer Dim factor As Double For i = 0 To n o(i) = i s(i) = Math.Abs(a(i, 0)) For j = 2 To n If Math.Abs(a(i, j)) > s(i) Then s(i) = Math.Abs(a(i, j)) End If Next j Next i For k = 0 To n - 1 Call Pivot(a, o, s, n, k) If Math.Abs(a(o(k), k) / s(o(k))) < tol Then er = -1 Exit For End If For i = k + 1 To n factor = a(o(i), k) / a(o(k), k) a(o(i), k) = factor For j = k + 1 To n a(o(i), j) = a(o(i), j) - factor * a(o(k), j) Next j Next i Next k If (Math.Abs(a(o(k), k) / s(o(k))) < tol) Then er = -1 End Sub Sub Pivot(a(,) As Double, o() As Double, s() As Double, n As int32, k As int32) Dim ii As Integer Dim p As Integer Dim big As Double Dim dummy As Double p = k big = Math.Abs(a(o(k), k) / s(o(k))) For ii = k + 1 To n dummy = Math.Abs(a(o(ii), k) / s(o(ii))) If dummy > big Then big = dummy p = ii End If Next ii dummy = o(p) o(p) = o(k) o(k) = dummy End Sub Sub Substitute(a(,) As Double, o() As Double, n As int32, b() As Double, x() As Double) Dim k As Integer Dim i As Integer Dim j As Integer Dim sum As Double Dim factor As Double For k = 0 To n - 1 For i = k + 1 To n factor = a(o(i), k) b(o(i)) = b(o(i)) - factor * b(o(k)) Next i Next k x(n) = b(o(n)) / a(o(n), n) For i = n - 1 To 0 Step -1 sum = 0 For j = i + 1 To n sum = sum + a(o(i), j) * x(j) Next j x(i) = (b(o(i)) - sum) / a(o(i), i) Next i End Sub

Consistently with the hipothesis that numerical methods for structural can be classifed into Time integrators, Matter Integrators, and Constraint Integrators, I have so far distinguished between LM (Lagrange Multipliers) and FEM (Finite Element Method), remaining for a near future the implemementation of a Newmark-Beta component.
The code provided above simply serves to calculate statics of structures, as the time is yet to be integrated.

I have made the convenient comparisons with a SAP2000 model to verify the exactness of the implementation, with very satisfactory results.

Well...more to come in future editions!