Stochastic Simulation and Lagrangian Dynamics: 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!

Stochastic Simulation and Lagrangian Dynamics

domingo, 25 de septiembre de 2011

First FEM + Lagrange multipliers program in Grasshopper

No hay comentarios:

Publicar un comentario

Statcounter

Where are we visited from?