Stochastic Simulation and Lagrangian Dynamics

A brief summary on the stiffness method for structural engineering

2012-01-31T14:02:00.000-08:00

Well, another year has gone by, and it has been a while since the last post. Life in Barcelona is being quite busy, so I am not dealing very well with my updates...This doesn't mean I am not researching, but apparently the topics I am getting to are deeper and deeper, so conclusions are getting harder to write...

The following is a summary of what I found here: www.civil.iitb.ac.in/~sghosh/CE317/history-stiffness-method.pdf

It is authored by the famous Zienkiewicz, so it must be a warranty of rigor. Out of it I have written my own summary, on pure chronological order, that goes as follows:

1862: Glebsch, Alfred wrote an algorithm for 3D trusses where no moments were involved. The models one would get out of it is as putting lots of springs connected by pins.
1883: Saint Venant, Andre wrote an article with comments on Glebsch's methodology.
1880: Manderla, Heinrich proposed an iterative method for the solution of the set of equations involved in pin-jointed trusses
1892: Mohr, Otto introduced the notion of "secondary stresses". He was referring to the moments. They were called secondary then as in steel and iron trusses the main ones were axial.
1910: Richardson provided the study of the stresses in a dam by means onf the Airy stress function.
1914: Bendixsen, Axel published a method similar to the slope deflection for structures with internal hyperstaticity.
1915: Wilson&Money proposed the slope deflection method, apparently unaware of Bendixsen's.
1922: Ĉaliŝev gave an iterative solution of successive approximations for frames without side-sway.
1930: Cross, Hardy published his famous method for redistributing unbalanced moments according to the internal stiffnesses of the composing bars. This provided the engineer with the necessary joint rotations.
1932: Grinter wrote a method for multi-story frames with side-sway called "method of successive corrections"
1935: Southwell tackled the problem of the constraints in a similar manner giving his famous method of "systematic relaxation of constraints", commonly known as dynamic relaxation.
1941: Courant, Richard explained his "variational methods for the solution of problems of equilibrium and vibrations". As he was a mathematician, the practical application of his work results still of some difficulty.
1943: McHenry in U.S. and Hrenikoff in Canada proposed a lattice analogy for the solution of stress problems so that continuum could be represented as a set of beams.
1944: Kron, Gabriel provided a complete algorithm in matrix form for 3D frames. In the 1930s it became fashionable to show the whole process of solving the linear equations in this systematic manner.
1954: Argyris wrote his "Energy theorems ans structural analysis", a comprehensive presentation of both the force and the displacement method for complicated airplane structures.
1967: Zienkiewicz presented his summary of the "Finite Element Method in structural Mechanics".

Tetrahedra discretization of a beam

2011-11-26T09:26:00.001-08:00

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

2011-10-31T14:42:00.000-07:00

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!

2011-10-09T09:44:00.000-07:00

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

2011-09-25T06:30:00.000-07:00

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!

Programming FEM with Grasshopper(R) for Rhino(R)

2011-07-17T09:04:00.000-07:00

In the past year a powerful design tool called Rhinoceros has come into my sphere of activity.
After several attempts of programming the necessary code to experiment with numerical methods using OGRE and C++, I have come to the conclusion that I need some more productive tools. The problem with C++ is that is good for final products and optimized code, but the available libraries for my purposes are normally too specific and complicated, and even simple debugging processes become a nightmare of libraries and strange variables.

Rhinoceros has an IDE plugin called Grasshopper that allows for visually programming with an object-oriented approach, which is very close to what we are intending now.
The next step in the research is to study the differences that appear when combining different methods in order to assess an optimum path for different applications. When studied with object-oriented mind, the numerical methods can easily be composed into GH components that afterwards we just need to link.

Arrangement of possible components with different integration concepts

One of my main concerns when developing with proprietary non-opensource IDEs and codes comes with the fact that there is always the chance of having to pay fees at some point, not to mention that they are opaque and not always supported by a community willing to help when there is trouble.
However, the case of Grasshopper seems quite strange as there is a strong and dynamic community of researchers and developers also often supported by excellent professionals (at least as far as I have researched, we'll see..).

Well, lets hope I have made the right decision and GH boosts my research results in the near future.

WHAT DOES SHAPE FUNCTION MEAN IN FINITE ELEMENT FORMULATION?

2011-07-17T08:05:00.000-07:00

Finite Element Analysis is a mathematical tool very extended among engineers. However, after more than a year researching on the topic of computer simulation, where FEA plays such an important role, I haven't yet found a satisfactory explanation on how they really really work...
Hopefully by means of this post I get to clarify some of the topics, trying to remove excessive algebraic and mathematical verbosity that is so annoyingly present everywhere.

The main background of FEM is that of structural engineering in the 60s. Engineers are very practical people, so initially they devised a system which allowed them to set algebraic equations where the relations between different points or "nodes" of their structures could be set.
These algebraic relations were further proven to be of many different types, not only structural, so also thermal relations could be formulated and many others: all one needed to analyze was a proper discretizacion of the space in the form of a mesh with nodes to relate to each other.

For our particular case of structural engineering (the one that matters for my PhD thesis), I have tried to illustrate the procedure in its three main steps, so the main ideas come up in a graphic taste:

STEP 1: Discretization
In this basic yet classic example, I have chosen to divide my sample structure in only 5 elements where nodes are clearly identifiable in the meetings between beams. There are four nodes.

Degrees of Freedom

Each node will have 6DOF (Degrees Of Freedom): three for linear displacement on each axis (X,Y,Z) and three for rotational around each axis (X,Y,Z), because we are working on 3D. Many examples available provide the more "simple" 2D situation, but in my opinion this only complicates things further.

A sample structure (mesh) and its topology table

Once the nodes are located and there is a network of how they relate to each other, we can consider we have a mesh. In our example, the correlation is depicted in the table: The elements serve to "link" nodes to each other. The table establishes a "topology" for the nodes.

STEP 2: Element characterization and shape functions
In this step is where FEM formulation and literature get really really awkward and nasty. In fact, this is the core of everything and where FEM differenciates from other ways of solving PDEs.

Different types of Finite Elements

In reality, and despite its mathematical complexity (also unneccesary to be explained so much in detail in my opinion), what we are looking for is a way of characterizing the material properties of the element. For such purpose, the method requires that the behavior of those links among nodes obeys some formula. This formula is the actual Shape Function. In fact, the shape function can be any mathematical formula that helps us to interpolate what happens wherever there are no points to define the mesh. This "ghost" entity that appears between nodes is in fact the Finite Element. In practical terms, as engineers we are more interested in the implementation of the FEM, not so much on its formulation, so what is important to understand is that for different shape functions we obtain different element matrices.

2DOF Beam element matrix

3DOF Beam Element matrix

3DOF Timoshenko Beam element matrix

6DOF Timoshenko Beam element matrix

Depending on the chosen formulation we have different degrees of interpolation and hence presumably higher or lower degrees of precision. Also depending on the chosen formulation we might have different ways of locating and relating our nodes to each other.

For our example we can choose any of the formulations provided in literature (above are the most common used in structural engineering). It is important to note the internal structure of these element matrices, which are symmetrical and clearly divided into parts each corresponding to the nodes that reside on the element's boundaries (2 nodes in the case of a beam - 4 quadrants in the matrix).

STEP 3: Matrix assembly and solution

Because the relations between nodes need to be accomplished all at the same time, we have to set all the equations in such a manner that they compose an algebraic system of equations. The matrix equation we want to solve (at least in statics) is as follows:

[F]= [Kg]·[u]

Where [F] is the vector of applied external forces, [Kg] is the system's global stiffness matrix and [u] is the vector of displacements resulting from the application of the forces. The size of [Kg], [F] and [u] is that of the number of DOFs times the number of nodes, being [Kg] squared and [F] and [u] unidimensional.

In order to generate the Force vector, all we need to do is collect the applied forces (linear and moment forces) and sort them according to the index of the node they are applied to.

For the global stiffness matrix, it is necessary a bit more laborious procedure by means of which we iterate throughout each element's particular stiffness matrix. Out of each one of those, we get only the part that corresponds to the position of the node we are storing in the matrix, and add it to the possible concurrent data that comes from other elements. Warning: before entering in the global stiffness matrix, we must convert local coordinates to global coordinates!

To get the displacement vector, it is needed to first enforce the constraints and then solve the resulting algebraic system of equations. To do this there are two classical approaches: Penalty Method and Lagrange Multipliers method, but we wont enter into this here...

Global Matrix Assembly

Afterwards, all we need to do is use any available algebraic equation solver (LU decomposition is one of the most extended), and obtain the solution of the system.

Impulses, FEM and Verlet

2011-05-09T09:20:00.000-07:00

After some time struggling with formulations of many diverse categories, I have eventually managed to implement a computer program where the scheme explained in the november entry (A proposal on a course on Real-Time Structural Dynamics) becomes clear and useful.The main point here is to explain how, by organizing the numerical methods according to what they really discretize, it is possible to grasp and understand such an entangled area of knowledge.One of the main difficulties that I am encountering through this research period is not just that of the thick mathematical language employed but mainly that of mixed concepts when authors come to explain each of the endless methods available in literature.The following is an step by step procedure where time, constraints and material properties are respectively integrated by means of Verlet, Impulse formulation and Finite Elements numerical methods.The software framework for development and visualization is that presented in April 2010 (one year already uuf!): First Ogre+Verlet+Gauss-Seidel simulation.On top of it I have made the necessary modifications and also the adaptation from Linux to Windows so it is easier to make further comparisons with commercial codes.The first that is needed is a set of nodes with their coordinates and their masses, along with the list of constraints that topologically relate one to another and to the environment:

Time integration: Verlet ODE integration method is known to have certain stability but not a big deal of accuracy. However, is fairly easy to implement. An exhaustive explanation on how it works can be found here: http://en.wikipedia.org/wiki/Verlet_integration. In our implementation, it is divided into functions: AccumulateAccelerations and Verlet.
In AccumulateAccelerations we iterate through each mass and add up earthquake readings from a file (per millisecond), gravity acceleration (980 cm/s2), and the internal forces caused by deformation divided by the weight of the mass.
The Verlet procedure, simply updates the position of each particle according to the formula:

x(t)=x(t-1) + v(t)·dt + a(t)·dt^2

Constraint integration: The only constraints included in the current implementation are those of the distance between two particles. By means of the concept of impulse, once two particles move close or apart by the effect of the accelerations, a corrective force is applied on each one of them. This force is applied in a very short lapse of time, hence can be regarded as an impulse. A more detailed explanation can be found in the work by Jan Bender (www.impulse-based.de). The algebraic system of equations that arises so as to satisfy all of the constraints simultaneously is solved by means of a Gauss-Seidel iterations.

Matter integration: For the consideration of matter properties the program uses the Finite Element Method. A quadratic shape function is used according to the implementation taught in http://www.colorado.edu/engineering/cas/courses.d/IFEM.d/. Chapters 20 and 21 explain clearly how to infer the values for the case of an elastic beam.
This provides our system with a rigidity matrix that is easily attached to each of the distance constraints. As the displacements have already been obtained previously, the solution of the force vector is straightforward by means of a simple rigidity matrix-displacements vector operation:

[f] = [K] · [u]

This force vector stores the effect of the deformation into the constraint and is later used to obtain the particle acceleration in the next timestep.

CANE STRUCTURES WITH FEM

2011-03-19T15:13:00.000-07:00

Currently Barcelona is holding a lot of activity in cane structure promoted by krfr collective (www.krfr.org). In a recent conversation with my colleague Oriol Palou (www.sustenta.eu) we discussed about the material properties of cane that are being studied within the School of Building Construction of Barcelona (EPSEB).

This led to a first experiment with FEM which video I show here. This simulates one load test made in EPSEB laboratories, where deflections in the middle of the arch under a 200 kg load were around 30 cm, with purely elastic behavior.
Obviously our model needs some refinement, as I had to load it with 2000 kg to make it deform that much, but the results are somehow encouraging

ON LAGRANGE'S UBIQUITI IN MODERN PHYSICS

2011-02-04T05:41:00.000-08:00

This post is just to point out a quite dazzling phenomenon I have been encountering during this year of research: Joseph Louis Lagrange.

It happens to be that this Italian gentleman revolutionized the world of Physics some two hundred years ago (see here a beautiful explanation on how), in such a manner that now is nearly impossible not to encounter his surname nearly everywhere when trying to understand them.
The following is a short outline of the mathematical/physical concepts including Lagrange (a larger version can be found in http://en.wikipedia.org/wiki/List_of_topics_named_after_Joseph_Louis_Lagrange):


Lagrange multipliers

Lagrange multipliers: These are mathematical artifacts for the solution of optimization problems (http://en.wikipedia.org/wiki/Lagrange_multiplier)
Euler-Lagrange equation: The Euler–Lagrange equation was developed in the 1750s by Euler and Lagrange in connection with their studies of the tautochrone problem. This is the problem of determining a curve on which a weighted particle will fall to a fixed point in a fixed amount of time, independent of the starting point. Lagrange solved this problem in 1755 and sent the solution to Euler. The two further developed Lagrange's method and applied it to mechanics, which led to the formulation of Lagrangian mechanics.
Euler-Lagrange equation
Their correspondence ultimately led to the calculus of variations, a term coined by Euler himself in 1766.In classical mechanics, it is equivalent to Newton's laws of motion, but it has the advantage that it takes the same form in any system of generalized coordinates, and it is better suited to generalizations. Is also related to optimization according to variational principles (http://en.wikipedia.org/wiki/Euler%E2%80%93Lagrange_equation).
Lagrangian function: The concept of a Lagrangian was originally introduced in a reformulation of classical mechanics by Irish mathematician William Rowan Hamilton known as Lagrangian mechanics.
Lagrangian function
In classical mechanics, the Lagrangian is defined as the kinetic energy, $T$ , of the system minus its potential energy.
Green-Lagrangian tensor: In continuum mechanics, the finite strain theory also called large strain theory, or large deformation theory, deals with deformations in which both rotations and strains are arbitrarily large. This means to invalidate the assumptions inherent in infinitesimal strain theory. In this case, the undeformed and deformed configurations of the continuum are significantly different and a clear distinction has to be made between them. This is commonly the case with elastomers, plastically-deforming materials and other fluids and biological soft tissue. The concept of strain is used to evaluate how much a given displacement differs locally from a rigid body displacement . One of such strains for large deformations is the Lagrangian finite strain tensor, also called the Green-Lagrangian strain tensor or Green - St-Venant strain tensor.
Lagrange description of motion. In continuum mechanics the Lagrangian specification of the flow field is a way of looking at fluid motion where the observer follows an individual fluid parcel as it moves through space and time. Plotting the position of an individual parcel through time gives the pathline of the parcel. This can be visualized as sitting in a boat and drifting down a river. The Eulerian specification of the flow field is a way of looking at fluid motion that focuses on specific locations in the space through which the fluid flows as time passes. This can be visualized by sitting on the bank of a river and watching the water pass the fixed location. The Lagrangian approach is also associated to particle based formulation, whereas the Eulerian is referred to as grid based formulations.
Eulerian description
Lagrangian description

And this is just a small sample. However, it fully justifies the name of the thesis and this blog's, as our research departs from Lagrange's work to try to find engineering solutions.
It must be said that Lagrange's prolificacy results somehow dazing and stunning, as there are so many fields he got involved in, and none of them of trivial nature.
I hope this quick outline serves others to find a way through all this tangled knowledge.

A proposal on a course on Real-Time Structural Dynamics

2010-11-22T14:08:00.000-08:00

The following is a tree with the key concepts that, although already at reach of any graduated structural engineer, need to be tied together, maybe in a course, in order to achieve a proper scope on how to simulate real-time structural dynamics:

Main time integration methods (ODE), their limitations (drawbacks), advantages, motivation, references, illustrations, year, all related to the three disciplines applied physics, maths and applied computing:

1st order

Euler
Backward Euler
Semi-implicit Euler
Exponential Euler

2nd order

Verlet
Velocity Verlet
Midpoint Method
Heun's
Newmark-beta
Leapfrog

Higher order

Runge-Kutta
Linear multistep

Main constraint/collision (DAE) integration methods, their limitations (drawbacks), advantages, motivation, references, illustrations, year, all related to the three disciplines applied physics, maths and applied computing.

Coordinate Partitioning
Constraint Orthogonalization
Udwadia-Kalaba

Main matter/continuum (PDE) integration methods, their limitations (drawbacks), advantages, motivation, references, illustrations, year, all related to the three disciplines applied physics, maths and applied computing.

MESH BASED METHODS

Finite Element
Finite Differences
Finite Volume
Boundary Element
Mass-spring systems

MESH FREE METHODS

SPH
Diffuse Element Method
Partition of Unity
Moving Least Square
Reproducing Kernel Method

A graphical approach to ths subject, not so tied to complex formalisms and formulation, would be a real help to attract more researchers into this fascinating discipline.
It also would reinforce the interest on the very differential equations, as these are a very abstract concept explained and taught on a very abstract basis. This makes them first candidate to be either forgotten or banned into the minds of students.

I will try to develop these subjects with more care in future issues...

Se vidimo!

Physics Engines Benchmarking

2010-10-26T06:49:00.000-07:00

The target was to devise in some manner the way different computationalphysics simulation models/engines/environments performed against acanonical one.

For such purpose, I have reviewed the paper "Beam Benchmark Problems For Validation Of Flexible Multibody Dynamic Codes" by A.L. Schwab and J.P. Meijgaard.
The following is an outline of the proposed in this paper:

1.- Introduction

The paper presents some basic problems for which analytic solution is known

2.- Beam benchmark problems

Tests to be performed

Static analyses in small displacements for the validation of the correct formulation of elastic forces
Static analyses for large displacements and rotations, on straight and curved beams
Bucklingtests in normal, lateral and torsional directions to check the waygeometric stiffness due to prestress is taken into account
Eigenfrequency analyses for the validation of the combination of elastic forces and distributed inertia forces
Mesh refinement tests for all the above to prove convergence of the results

Underlying model: Timoshenko beam with large displacements and large rotations

Finite Element Method Beam Element
Shear flexible based on the elastic line concept
Slender beam, cross section doubly symmetric
Large rotation and displacements, but small deformations
Isotropic and linearly elastic
BEAM model: standard strain-displacement relations
BEAMNL model: additional quadratic terms included in the strain-displacement for better performance in the pre-stress cases

Interestingly enough, the analytical solution isprovided for each of the tests, and the proposed model is simple enoughso as not having too much trouble in introducing it in a fewcharacteristic ready made softwares currently available.

The implementation of such a benchmark over some of the most popular engines revised lately would surely make a nice contribution.

Integration Overview

2010-10-14T04:49:00.000-07:00

In order to get some scope I have prepared a new diagram where the different integration fields of any multibody physics engine can be fit.

Normally we have to integrate time, via ODEs, for which any of the available schemes can be chosen, but then also the Differential Algebraic Equations for constraints (related to collisions) and the Partial Differential Equations related to the continuum have to be solved.

Of course, not every engine implements the continuum part (limiting to rigid solid) and some matter integration schemes already consider the very collisions so the constraint integration is explicitly sorted...

The associated disciplines where each concept fits are represented by the horizontal blocks.

Discipline Overview

2010-09-22T04:45:00.000-07:00

Holidays are over already. I am now living in Slovenia, Ljubljana, and developing the thesis in the Faculty of Architecture (lovely place, by the way...).

Last June I presented my Thesis Proposal to the jury of the Computational and Applied Physics Department in the UPC and...it was accepted!!

After the long pause I have had to prepare an overview of all the studied subjects, just to get some scope on what is coming and where we are now. The following diagrams give some idea on how it is:

And with a bit more of detail:

These diagrams derive from the need of disentangling disciplines that more often than seldom get mixed in the literature and make it fairly hard to follow.
The literature I am talking about is that of the computational simulation of physics, where numerical methods meet with very different origins.

SUSTAINABLE DYNAMIC STRUCTURES WITH CANE

2010-07-27T08:45:00.000-07:00

The last month I have been quite busy helping out raise the structure above.
It is located near Idanha Velha, in the premises for the International Boom Festival in Portugal.

It is one of the most clear examples of dynamic structures, where ordinary laws of statics don't apply, and one of my first fields of intense research.

The aim is to reach a comprehensive computational model that allows the designer (Jonathan Cory Wright - www.canyaviva.com) to experiment with his designs prior to having to actually raise his buildings!

An Overview on Constraint Enforced Formulations of Variational Dynamics

2010-04-20T14:55:00.000-07:00

A couple of weeks ago a very interesting paper entered my hard drive: Review of classical approaches for constraint enforcement in multibody systems.
It is a fairly clarifying overview on solid grounds of certain methodologies to solve multibody systems dynamics.
These methodologies differ from the ones I have encountered in the spreaded dynamics engines in their apparently more robust and simple formulation.
Here is a synthesis of this paper:

Multibody systems present two distinguishable features:
Bodies undergo finite relative rotations, which introduce nonlinearities
Bodies are connected by mechanical joints that impose restrictions, which mean a set of governing equations that combine differential and algebraic equations (ODEs and DAEs, respectively).
Lagrange's equation of the first kind has the following aspect:
Where
M=M(q,t) is the mass matrix
q is the generalized coordinates vector
B is the constraint matrix
λ is the array of Lagrange multipliers
F is the dynamic externally applied forces
If all the constraints are holonomic (velocity independent), B is called the Jacobian matrix, and the generalized coordinates, q, are linked by m algebraic constraints.
Lagrange's equations of the first kind form a set of (m+n) Differential Algebraic Equations.
The approach of the following methods is to use algebraic procedures in order to eliminate Lagrange's multipliers and then obtain a set of ODEs.
Maggi's formulation: implies the creation of a vector containing the so called kinematic parameters, generalized speeds or independent quasi-velocities by the analyst in order to obtain a Γ matrix. This matrix spans the null space of the constraint matrix and allows for the elimination of the Lagrange multipliers.
Index-1 formulation: requires that initial condition of the problem be subjected to the constraint conditions. Then, it is possible to obtain a system of 2nd order ODEs that is solvable by rearranging the previous equation, extracting the Lagrange multipliers and hence obtaining:
Null space formulation: this method solves the system of second order ODEs by premultiplying the first part of the equation by the transposed null space matrix thus eliminating the Lagrange multipliers.
Udwadia-Kalaba formulation: this method represents a more compact and general form of solving the DAEs by means of the Moore-Penrose generalized inverse. It is based on Gauss' Principle of Minimum Constraint, which establishes that the explicit equations of motion be expressed as the solution of a quadratic minimization problem subjected to constraints, but at the acceleration level.

All these formulations transform the (2n+m) first order DAEs into ODEs by eliminating Lagrange multipliers.
Maggi's formulation yields (2n-m) first order ODEs.
Index-1, null space and Udwadia-Kalaba form sets of (n) second order ODEs which could be alternatively recast into (2n) first order ODEs for the n generalized coordinates and the n generalized velocities.
The main advantage of these methodologies is not so much the reduction in the number of equations but rather in the change from DAEs to ODEs.
There is a warning on the constraint drift phenomenon, for which these method will be more affected (not so much in Maggi's formulation), and that would require constraint stabilization techniques.

My conclusion is that it seems the way to go, not only for the claims of more stable and quick numerical techniques available to get them working, but also for an apparently more clear approach in the theoretical field.
Particularly, I have done some research into the Udwadia-Kalaba formulation, and definitively is a very promising one.

ODE DYNAMICS ENGINE QUICK OVERVIEW

2010-04-15T15:26:00.000-07:00

After seeing my own simulator running there is a whole lot of things that I wish it had.
All this stuff is more or less complicated to achieve, but from my previous research I have found some Open Source engines from which I can learn a lot.

The best documented one so far is ODE (www.ode.org), and also the one I intend to merge with Blender.
So, I present a very quick and thorough overview on how this works, in order to have some scope.

The first thing to disentangle is how it gets working. For that, just looking at the user's manual one gets the following algorithm:

1.-Create a Dynamics world
This basically means to make an instance of the dxWorld struct.
2.-Create the bodies
Attach the bodies to the world (this means adding data to the dxBody* array of the dynamics world).
Set their properties (position and orientation of point of reference, its linear and angular velocities, the mass of the bodies and some other stuff of the like).
3.-Create the joints
Attach the joints to the world (by adding data to the dxJoints* array of the dynamics world)
Set their properties (depending of the selected type of joint, one has to provide different details).
4.-Manage collisions
Create a new collision world (just by making an instance to the dxSpace struct).
Create a joint group where collisions will be stored temporarily for every frame step.
5.-Loop
Apply forces to bodies.
Adjust joint parameters.
Call collision detection.
Create a contact joint for every detected collision point and add it to the collision joint group.
Take a simulation step.
Clear the collision joint group.
6.-Destroy the dynamics and the collision worlds (wow, that sounds evil...hehehe).

All of this seems fairly easy to do but then one has to get to know where and how things have been implemented. This is when things become a bit harder.
I am sure the file structure and the class definitions make perfect sense for the programmers of this engine.
Fortunately, they have been careful enough so as to comment everything in an understandable way, for which I feel deeply grateful.
After some digging, I have managed to restructure the ode/src folder into the following topics:

Accessories
Memory management
Math
Matrix handling
External applications
Collision
Core
Joints
Primitives
Solver

Which means one can tackle the engine in an ordered manner and find things when needed.

From a theoretical point of view, this engine presents the following features:

A Lagrange multiplier velocity based model from Trinkle and Stewart
A friction model from Baraff
A Danzig LCP solver

First Ogre+Verlet+Gauss-Seidel simulation

2010-04-13T20:14:00.000-07:00

Once the Ogre3D engine is ready to draw what we numerically compute, it was about time to start having some fun.

I found this excellent article from Thomas Jakobsen: http://www.gotoandplay.it/_articles/2005/08/advCharPhysics.php where a very neat and simple engine is implemented (thanks Mr. Jakobsen).

It is very well explained so adapting it to my little lab was not very hard.
Here I proudly present my very first simulation!

I am perfectly conscious that it lacks of a lot of things...namely collision detection, a proper stable, precise integrator, an optimized implementation...and a long etcetera...but it's just a baby!

Here is a zip file with the code (an extension from what I presented in the last post):

http://www.mnbvlabs.com/Thesis/VerletRelaxation.zip

Setting Ogre3D as Graphic Environment

2010-04-11T15:51:00.000-07:00

The last days I have been working on getting a Graphics Environment for my Toy Physics Engine.
Last December and January I had already fiddled around with excellent engines Ogre (www.ogre3d.org) and Bullet (http://bulletphysics.org/wordpress/), trying to put them to work together.I have a post from february where some simulations on youtube can be seen and some preliminary introductions are made.

However, now I am starting to understand things, and whenever I have the drive to test them always encounter the problem on how visualize them.

So, I came across with Blender (www.blender.org), another fantastic tool, and got surprised to find their Google Summer of Code mentorship for 2010 (http://wiki.blender.org/index.php/Dev:Ref/GSoC/2010/Info).
I have applied for it. Deadline was 9th April, which meant preparing everything in a hurry and some embarrasing trouble with their wiki page (I have already apologised in their bf-commiters mailing list...ooops).
Here is my proposal for such an important event: (http://wiki.blender.org/index.php/User:Ndujar).
Hopefully they find it interesting and contact me to complete it!

In the meantime, I have prepared a short and basic Ogre framework to begin doing some tests.
It is composed of a few files and classes integrating Ogre.

I generated it basically copying, cutting and pasting from the very good tutorials available at the Ogre web page (http://www.ogre3d.org/wiki/index.php/Ogre_Tutorials).

It goes as follows:

The first step is to intall and cofigure Ogre in your computer. This is the best tutorial I found: http://ubuntuforums.org/archive/index.php/t-1148570.html.

Then, using whatever IDE of your taste (I am on EasyEclipse), just get the following files to compile.

The core of the application is the BaseApp object, from which I will later extend further as needed. It is basically in charge of doing everything Ogre requires to be done to put things up in a rendering window.
BaseApp.h:

BaseApp.cpp:

#include "BaseApp.h"using namespace Ogre;using namespace std; ///This code is adapted from the Ogre tutorials webpage for my own purposes ///http://www.ogre3d.org/wiki/index.php/Basic_Tutorial_6#Creating_the_RenderSystem ///Main function, calling all the private settings to make the application run in proper order: ///The basic Ogre life cycle looks like this: ///   1. Create the Root object. ///   2. Define the resources that Ogre will use. ///   3. Choose and set up the RenderSystem (that is, DirectX, OpenGL, etc). ///   4. Create the RenderWindow (the window which Ogre resides in). ///   5. Initialize the resources that you are going to use. ///   6. Create a scene using those resources. ///   7. Set up any third party libraries and plugins. ///   8. Create any number of frame listeners. ///   9. Start the render loop.    void BaseApp::go()    {        createRoot();        defineResources();        loadResources();        setupRenderSystem();        createRenderWindow();        initializeResourceGroups();        setupScene();        createFrameListener();        createPhysicsFrameListener();        startRenderLoop();    }    ///Destructor releases the FrameListener and the Root object from memory    BaseApp::~BaseApp()    {        delete mListener;        delete mRoot;    } ///-createRoot sets the Ogre Root object from where all others derive.    ///The Root object is the core of the Ogre library,    ///and must be created before you can do almost anything with the engine.    void BaseApp::createRoot()    {        mRoot = new Root();    } ///-defineResources sets all the Ogre engine requirements by iterating the configuration file.    ///The next thing we have to do is define the resources that the application uses.    ///This includes the textures, models, scripts, and so on.    ///Everything is set up in the resources.cfg file.    void BaseApp::defineResources()    // Method which will define the source of resources (other than current folder)    {        // Load resource paths from config file        ConfigFile cf;        mResourcePath = "";        cf.load(mResourcePath + "resources.cfg");        // Go through all sections & settings in the file        ConfigFile::SectionIterator seci = cf.getSectionIterator();        String secName, typeName, archName;        while (seci.hasMoreElements())        {            secName = seci.peekNextKey();            ConfigFile::SettingsMultiMap *settings = seci.getNext();            ConfigFile::SettingsMultiMap::iterator i;            for (i = settings->begin(); i != settings->end(); ++i)            {                typeName = i->first;                archName = i->second;                ResourceGroupManager::getSingleton().addResourceLocation(archName, typeName, secName);            }        }    } ///-loadResources intializes the application ResourceGroupManager.    ///The function defineResources adds all of the resources from the config file,    ///but it only tells Ogre where they are.    ///Before you can use any of them you have to either initialize the group that you want to use,    ///or you must initialize them all.    void BaseApp::loadResources()    {       // Initialise, parse scripts etc       ResourceGroupManager::getSingleton().initialiseAllResourceGroups();    } ///-setupRenderSystem checks wether the config file Ogre.cfg has already set things up and calls the /// Ogre render system setup dialog in case it isn't.    ///In the first part of the if statement, we attempt to restore the config file.    ///If that function returns false it means that the file does not exist so we should show the config dialog,    ///which is the second portion of that if statement.    ///If that also returns false it means the user canceled out of config dialog    ///(meaning they want to exit the program).    void BaseApp::setupRenderSystem()    {        if (!mRoot->restoreConfig() && !mRoot->showConfigDialog())            throw Exception(52, "User canceled the config dialog!", "BaseApp::setupRenderSystem()");    } ///-createRenderWindow sets the Ogre root window where everything is going to be displayed    ///Now that we have chosen the RenderSystem, we need a window to render Ogre in.    ///There are actually a lot of options for how to do this, but we will really only need this to    ///keep thing simple.    void BaseApp::createRenderWindow()    {        mRoot->initialise(true, "Physics World");    }    ///The only thing left to do is to initialize the resources we are about to use.    ///In a very large game or application, we may have hundreds or even thousands of    ///resources that our game uses - everything from meshes to textures to scripts.    ///At any given time though, we probably will only be using a small subset of these resources.    ///To keep down memory requirements, we can load only the resources that our application is using.    ///We do this by dividing the resources into sections and only initializing them as we go.    void BaseApp::initializeResourceGroups()    {        TextureManager::getSingleton().setDefaultNumMipmaps(5);        ResourceGroupManager::getSingleton().initialiseAllResourceGroups();    }    ///There are at least three things that must be done before you start adding things to a scene:    ///creating the SceneManager, creating the Camera, and creating the Viewport.    ///I have added some lights and entities to give it some flesh.    ///Entities are instances of PhysicsEntity, which is the class where all physics properties    ///will be implemented    void BaseApp::setupScene()    {        mgr = mRoot->createSceneManager(ST_GENERIC, "Default SceneManager");        Viewport *vp;        PhysicsEntity *PhysEnt(0);        Light *light;        mgr->setAmbientLight( ColourValue( 0, 0, 0 ) );        mgr->setShadowTechnique( SHADOWTYPE_TEXTURE_MODULATIVE );        mCamera = mgr->createCamera("PlayerCam");        mCamera->setPosition(Vector3(300,300,200));        mCamera->lookAt(Vector3(0,0,0));        mCamera->setNearClipDistance(5);        vp = mRoot->getAutoCreatedWindow()->addViewport(mCamera);        vp->setBackgroundColour(ColourValue(0,0,0));        mCamera->setAspectRatio(Real(vp->getActualWidth()) / Real(vp->getActualHeight()));        light = mgr->createLight("Light1");        light->setType(Light::LT_POINT);        light->setPosition(0, 15, 25);        light->setDiffuseColour(.25, .25, 0);        light->setSpecularColour(.25, .25, 0);        light = mgr->createLight("Light3");        light->setType(Light::LT_DIRECTIONAL);        light->setPosition(10,10,10);        light->setCastShadows(true);        light->setDiffuseColour(.25, .25, 0);        light->setSpecularColour(.25, .25, 0);        light->setDirection(Vector3( 0, -1, 1 ));        PhysEnt=new PhysicsEntity(mgr,PhysicsEntity::cube);        PhysEnt->Geometry->scale(10,10,10);        PhysicsEntities.push_back(*PhysEnt);        PhysEnt=new PhysicsEntity(mgr,PhysicsEntity::plane);        PhysEnt->Geometry->scale(1000,1000,1000);        PhysicsEntities.push_back(*PhysEnt);    }    ///Every frame you must call the capture method on all Keyboard, Mouse, and Joystick objects you use.    ///Our BaseFrameListener class manages all this for us.    void BaseApp::createFrameListener()    {     mListener= new BaseFrameListener(mRoot->getAutoCreatedWindow(),mCamera);        mRoot->addFrameListener(mListener);    }    ///Every frame you must also update all physics entities    ///Our PhysicsFramelistener class manages all this for us.    void BaseApp::createPhysicsFrameListener()    {     PhListener= new PhysicsFrameListener(PhysicsEntities);        mRoot->addFrameListener(PhListener);    }    ///This will render the application until a FrameListener returns false.    void BaseApp::startRenderLoop()    {        mRoot->startRendering();    }

The BaseFrameListener is another very important piece of code that allows one to control the rendering flow.
In order to manage whatever happens on every rendering step, Ogre provides the FrameListener object.
With the installed package comes an ExampleFrameListener.h file, which overrides the basic Ogre FrameListener and allows for a fairly good control over the inputting devices such as keyboard and mouse.I have slightly customised it for my purposes.
It has been a bit tricky however, and has made my computer crash a few times.
Anyway, here it is, domesticated at last!

BaseFrameListener.h:

/*-----------------------------------------------------------------------------This source file is part of OGRE    (Object-oriented Graphics Rendering Engine)For the latest info, see http://www.ogre3d.org/Copyright (c) 2000-2006 Torus Knot Software LtdAlso see acknowledgements in Readme.htmlYou may use this sample code for anything you like, it is not covered by theLGPL like the rest of the engine.-----------------------------------------------------------------------------*//*-----------------------------------------------------------------------------Filename:    ExampleFrameListener.h->Adapted into BaseFrameListenerDescription: Defines an example frame listener which responds to frame events.This frame listener just moves a specified camera around based onkeyboard and mouse movements.Mouse:    FreelookW or Up:  ForwardS or Down:BackwardA:        Step leftD:        Step rightPgUp:     Move upwardsPgDown:   Move downwardsF:        Toggle frame rate stats on/offR:        Render modeT:        Cycle texture filtering                       Bilinear, Trilinear, Anisotropic(8)P:        Toggle on/off display of camera position / orientationImp Pant: Saves an incremental Screenshot-----------------------------------------------------------------------------*/#ifndef BASEFRAMELISTENER_H_#define BASEFRAMELISTENER_H_#include "Ogre.h"#include "OgreStringConverter.h"#include "OgreException.h"//Use this define to signify OIS will be used as a DLL//(so that dll import/export macros are in effect)#define OIS_DYNAMIC_LIB#include <OIS/OIS.h>using namespace Ogre;class BaseFrameListener: public FrameListener, public WindowEventListener{protected: Camera* mCamera; Vector3 mTranslateVector; Real mCurrentSpeed; RenderWindow* mWindow; bool mStatsOn; std::string mDebugText; unsigned int mNumScreenShots; float mMoveScale; float mSpeedLimit; Degree mRotScale; // just to stop toggles flipping too fast Real mTimeUntilNextToggle ; Radian mRotX, mRotY; TextureFilterOptions mFiltering; int mAniso; int mSceneDetailIndex ; Real mMoveSpeed; Degree mRotateSpeed; Overlay* mDebugOverlay; //OIS Input devices OIS::InputManager* mInputManager; OIS::Keyboard* mKeyboard; OIS::Mouse*    mMouse; virtual void updateStats(void) {  static String currFps = "Current FPS: ";  static String avgFps = "Average FPS: ";  static String bestFps = "Best FPS: ";  static String worstFps = "Worst FPS: ";  static String tris = "Triangle Count: ";  static String batches = "Batch Count: ";  // update stats when necessary  try {   OverlayElement* guiAvg = OverlayManager::getSingleton().getOverlayElement("Core/AverageFps");   OverlayElement* guiCurr = OverlayManager::getSingleton().getOverlayElement("Core/CurrFps");   OverlayElement* guiBest = OverlayManager::getSingleton().getOverlayElement("Core/BestFps");   OverlayElement* guiWorst = OverlayManager::getSingleton().getOverlayElement("Core/WorstFps");   const RenderTarget::FrameStats& stats = mWindow->getStatistics();   guiAvg->setCaption(avgFps + StringConverter::toString(stats.avgFPS));   guiCurr->setCaption(currFps + StringConverter::toString(stats.lastFPS));   guiBest->setCaption(bestFps + StringConverter::toString(stats.bestFPS)    +" "+StringConverter::toString(stats.bestFrameTime)+" ms");   guiWorst->setCaption(worstFps + StringConverter::toString(stats.worstFPS)    +" "+StringConverter::toString(stats.worstFrameTime)+" ms");   OverlayElement* guiTris = OverlayManager::getSingleton().getOverlayElement("Core/NumTris");   guiTris->setCaption(tris + StringConverter::toString(stats.triangleCount));   OverlayElement* guiBatches = OverlayManager::getSingleton().getOverlayElement("Core/NumBatches");   guiBatches->setCaption(batches + StringConverter::toString(stats.batchCount));   OverlayElement* guiDbg = OverlayManager::getSingleton().getOverlayElement("Core/DebugText");   guiDbg->setCaption(mDebugText);  }  catch(...) { /* ignore */ } }public: // Constructor takes a RenderWindow because it uses that to determine input context BaseFrameListener(RenderWindow* win, Camera* cam, bool bufferedKeys = false, bool bufferedMouse = false) :  mCamera(cam),  mTranslateVector(Vector3::ZERO),  mCurrentSpeed(0),  mWindow(win),  mStatsOn(true),  mNumScreenShots(0),  mMoveScale(0.0f),  mRotScale(0.0f),  mTimeUntilNextToggle(0),  mFiltering(TFO_BILINEAR),  mAniso(1),  mSceneDetailIndex(0),  mMoveSpeed(100),  mRotateSpeed(36),  mDebugOverlay(0),  mInputManager(0),  mKeyboard(0),  mMouse(0) {  mDebugOverlay = OverlayManager::getSingleton().getByName("Core/DebugOverlay");  LogManager::getSingletonPtr()->logMessage("*** Initializing OIS ***");  OIS::ParamList pl;  size_t windowHnd = 0;  std::ostringstream windowHndStr;  win->getCustomAttribute("WINDOW", &windowHnd);  windowHndStr << windowHnd;  pl.insert(std::make_pair(std::string("WINDOW"), windowHndStr.str()));  mInputManager = OIS::InputManager::createInputSystem( pl );  //Create all devices (We only catch joystick exceptions here, as, most people have Key/Mouse)  mKeyboard = static_cast<OIS::Keyboard*>(mInputManager->createInputObject( OIS::OISKeyboard, bufferedKeys ));  mMouse = static_cast<OIS::Mouse*>(mInputManager->createInputObject( OIS::OISMouse, bufferedMouse ));  //Set initial mouse clipping size  windowResized(mWindow);  showDebugOverlay(true);  //Register as a Window listener  WindowEventUtilities::addWindowEventListener(mWindow, this); } //Adjust mouse clipping area virtual void windowResized(RenderWindow* rw) {  unsigned int width, height, depth;  int left, top;  rw->getMetrics(width, height, depth, left, top);  const OIS::MouseState &ms = mMouse->getMouseState();  ms.width = width;  ms.height = height; } //Unattach OIS before window shutdown (very important under Linux) virtual void windowClosed(RenderWindow* rw) {  //Only close for window that created OIS (the main window in these demos)  if( rw == mWindow )  {   if( mInputManager )   {    mInputManager->destroyInputObject( mMouse );    mInputManager->destroyInputObject( mKeyboard );    OIS::InputManager::destroyInputSystem(mInputManager);    mInputManager = 0;   }  } } virtual ~BaseFrameListener() {  //Remove ourself as a Window listener  WindowEventUtilities::removeWindowEventListener(mWindow, this);  windowClosed(mWindow); } virtual bool processUnbufferedKeyInput(const FrameEvent& evt) {  if(mKeyboard->isKeyDown(OIS::KC_A))   mTranslateVector.x = -mMoveScale; // Move camera left  if(mKeyboard->isKeyDown(OIS::KC_D))   mTranslateVector.x = mMoveScale; // Move camera RIGHT  if(mKeyboard->isKeyDown(OIS::KC_UP) || mKeyboard->isKeyDown(OIS::KC_W) )   mTranslateVector.z = -mMoveScale; // Move camera forward  if(mKeyboard->isKeyDown(OIS::KC_DOWN) || mKeyboard->isKeyDown(OIS::KC_S) )   mTranslateVector.z = mMoveScale; // Move camera backward  if(mKeyboard->isKeyDown(OIS::KC_PGUP))   mTranslateVector.y = mMoveScale; // Move camera up  if(mKeyboard->isKeyDown(OIS::KC_PGDOWN))   mTranslateVector.y = -mMoveScale; // Move camera down  if(mKeyboard->isKeyDown(OIS::KC_RIGHT))   mCamera->yaw(-mRotScale);  if(mKeyboard->isKeyDown(OIS::KC_LEFT))   mCamera->yaw(mRotScale);  if( mKeyboard->isKeyDown(OIS::KC_ESCAPE) || mKeyboard->isKeyDown(OIS::KC_Q) )   return false;        if( mKeyboard->isKeyDown(OIS::KC_F) && mTimeUntilNextToggle <= 0 )  {   mStatsOn = !mStatsOn;   showDebugOverlay(mStatsOn);   mTimeUntilNextToggle = 1;  }  if( mKeyboard->isKeyDown(OIS::KC_T) && mTimeUntilNextToggle <= 0 )  {   switch(mFiltering)   {   case TFO_BILINEAR:    mFiltering = TFO_TRILINEAR;    mAniso = 1;    break;   case TFO_TRILINEAR:    mFiltering = TFO_ANISOTROPIC;    mAniso = 8;    break;   case TFO_ANISOTROPIC:    mFiltering = TFO_BILINEAR;    mAniso = 1;    break;   default: break;   }   MaterialManager::getSingleton().setDefaultTextureFiltering(mFiltering);   MaterialManager::getSingleton().setDefaultAnisotropy(mAniso);   showDebugOverlay(mStatsOn);   mTimeUntilNextToggle = 1;  }  if(mKeyboard->isKeyDown(OIS::KC_SYSRQ) && mTimeUntilNextToggle <= 0)  {   std::ostringstream ss;   ss << "screenshot_" << ++mNumScreenShots << ".png";   mWindow->writeContentsToFile(ss.str());   mTimeUntilNextToggle = 0.5;   mDebugText = "Saved: " + ss.str();  }  if(mKeyboard->isKeyDown(OIS::KC_R) && mTimeUntilNextToggle <=0)  {   mSceneDetailIndex = (mSceneDetailIndex+1)%3 ;   switch(mSceneDetailIndex) {    case 0 : mCamera->setPolygonMode(PM_SOLID); break;    case 1 : mCamera->setPolygonMode(PM_WIREFRAME); break;    case 2 : mCamera->setPolygonMode(PM_POINTS); break;   }   mTimeUntilNextToggle = 0.5;  }  static bool displayCameraDetails = false;  if(mKeyboard->isKeyDown(OIS::KC_P) && mTimeUntilNextToggle <= 0)  {   displayCameraDetails = !displayCameraDetails;   mTimeUntilNextToggle = 0.5;   if (!displayCameraDetails)    mDebugText = "";  }  // Print camera details  if(displayCameraDetails)   mDebugText = "P: " + StringConverter::toString(mCamera->getDerivedPosition()) +       " " + "O: " + StringConverter::toString(mCamera->getDerivedOrientation());  // Return true to continue rendering  return true; } virtual bool processUnbufferedMouseInput(const FrameEvent& evt) {  // Rotation factors, may not be used if the second mouse button is pressed  // 2nd mouse button - slide, otherwise rotate  const OIS::MouseState &ms = mMouse->getMouseState();  if( ms.buttonDown( OIS::MB_Right ) )  {   mTranslateVector.x += ms.X.rel * 0.13;   mTranslateVector.y -= ms.Y.rel * 0.13;  }  else  {   mRotX = Degree(-ms.X.rel * 0.13);   mRotY = Degree(-ms.Y.rel * 0.13);  }  return true; } virtual void moveCamera() {  // Make all the changes to the camera  // Note that YAW direction is around a fixed axis (freelook style) rather than a natural YAW  //(e.g. airplane)  mCamera->yaw(mRotX);  mCamera->pitch(mRotY);  mCamera->moveRelative(mTranslateVector); } virtual void showDebugOverlay(bool show) {  if (mDebugOverlay)  {   if (show)    mDebugOverlay->show();   else    mDebugOverlay->hide();  } } // Override frameRenderingQueued event to process that (don't care about frameEnded) bool frameRenderingQueued(const FrameEvent& evt) {  if(mWindow->isClosed()) return false;  mSpeedLimit = mMoveScale * evt.timeSinceLastFrame;  //Need to capture/update each device  mKeyboard->capture();  mMouse->capture();     Ogre::Vector3 lastMotion = mTranslateVector;  //Check if one of the devices is not buffered  if(!mKeyboard->buffered()|| !mMouse->buffered() )  {   // one of the input modes is immediate, so setup what is needed for immediate movement   if (mTimeUntilNextToggle >= 0)    mTimeUntilNextToggle -= evt.timeSinceLastFrame;   // Move about 100 units per second   mMoveScale = mMoveSpeed * evt.timeSinceLastFrame;   // Take about 10 seconds for full rotation   mRotScale = mRotateSpeed * evt.timeSinceLastFrame;   mRotX = 0;   mRotY = 0;   mTranslateVector = Ogre::Vector3::ZERO;  }  //Check to see which device is not buffered, and handle it  if( !mKeyboard->buffered() )   if( processUnbufferedKeyInput(evt) == false )    return false;  if( !mMouse->buffered() )   if( processUnbufferedMouseInput(evt) == false )    return false;  // ramp up / ramp down speed     if (mTranslateVector == Ogre::Vector3::ZERO)  {   // decay (one third speed)   mCurrentSpeed -= evt.timeSinceLastFrame * 0.3;   mTranslateVector = lastMotion;  }  else  {   // ramp up   mCurrentSpeed += evt.timeSinceLastFrame;  }  // Limit motion speed  if (mCurrentSpeed > 1.0)   mCurrentSpeed = 1.0;  if (mCurrentSpeed < 0.0)   mCurrentSpeed = 0.0;  mTranslateVector *= mCurrentSpeed;  if(!mKeyboard->buffered() || !mMouse->buffered())   moveCamera();  return true; } bool frameEnded(const FrameEvent& evt) {  updateStats();  return true; }};#endif /*BASEFRAMELISTENER_H_*/

Once the beast is under control, the rest is a piece of cake...hehehe...
You can make as many instances of a FrameListener as you need in an Ogre application.
I have chosen to create a separate framelistener for managing whatever goes under the physics, to keep things up neat and tidy.
It is the PhysicsFrameListener class:

PhysicsFrameListener.h:

/*#ifndef PHYSICSFRAMELISTENER_H_#define PHYSICSFRAMELISTENER_H_#include "PhysicsEntity.h"#include <iostream>#include <math.h>#include <vector>#include <time.h>#include <ctime>#include <stdio.h>#include <unistd.h>#include <stdlib.h>#include <Ogre.h>using namespace Ogre;///This class is in charge to manage the physics before every///render step on the Ogre window.class PhysicsFrameListener : public FrameListener{public:    ///The Ents vector contains the entities with physical properties attached to them    ///In this example, no physics are used yet. Just an X-axis displacement is applied    ///to perceive motion.    PhysicsFrameListener(std::vector<PhysicsEntity> Ents):     Entities(Ents)    {    }    ///frameStarted is an overriden event from the Ogre FrameListener.    ///Note that positions are updated before sending the information    ///to the root object in the go loop.    bool frameStarted(const FrameEvent& evt)    {     UpdatePositions();     return true;    }    ///Update positions is in charge to reallocate the entities to their future locations    void UpdatePositions()    {//     PhysicsEntity::EntityType TypeEnt;     for (unsigned int i=0; i<Entities.size(); i++)     {      if (Entities[i].GeomType==PhysicsEntity::cube)       Entities[i].Geometry->translate(0.1,0,0,Ogre::Node::TS_WORLD);     }    }private:    std::vector<PhysicsEntity> Entities;};#endif /*PHYSICSFRAMELISTENER_H_*/

The elementary brick on which everything else will be constructed is likely to be called PhysicsEntity.
So I have defined my PhysicsEntity Class.
In this code is still at its very bones, just presenting a couple of functions to make things visible.
Hopefully in the near future some more flesh will get added on it!

PhysicsEntity.h:

PhysicsEntity.cpp:

/*#include "PhysicsEntity.h"PhysicsEntity::PhysicsEntity(SceneManager *SceneMgr,EntityType ShapeType){ Geometry = SceneMgr->getRootSceneNode()->createChildSceneNode(); Geometry->setPosition(0,5,0); GeomType=ShapeType; switch (GeomType) { case cube:        Geometry->attachObject(createCubeMesh("Cube","myMaterial"));  break; case plane:        Geometry->attachObject(createPlane("Plane","myMaterial"));  break; };}PhysicsEntity::~PhysicsEntity(){}ManualObject* PhysicsEntity::createCubeMesh(Ogre::String name, Ogre::String matName){ ManualObject* cube = new ManualObject(name); cube->begin(matName); cube->position(0.5,-0.5,1.0);cube->normal(0.408248,-0.816497,0.408248);cube->textureCoord(1,0); cube->position(-0.5,-0.5,0.0);cube->normal(-0.408248,-0.816497,-0.408248);cube->textureCoord(0,1); cube->position(0.5,-0.5,0.0);cube->normal(0.666667,-0.333333,-0.666667);cube->textureCoord(1,1); cube->position(-0.5,-0.5,1.0);cube->normal(-0.666667,-0.333333,0.666667);cube->textureCoord(0,0); cube->position(0.5,0.5,1.0);cube->normal(0.666667,0.333333,0.666667);cube->textureCoord(1,0); cube->position(-0.5,-0.5,1.0);cube->normal(-0.666667,-0.333333,0.666667);cube->textureCoord(0,1); cube->position(0.5,-0.5,1.0);cube->normal(0.408248,-0.816497,0.408248);cube->textureCoord(1,1); cube->position(-0.5,0.5,1.0);cube->normal(-0.408248,0.816497,0.408248);cube->textureCoord(0,0); cube->position(-0.5,0.5,0.0);cube->normal(-0.666667,0.333333,-0.666667);cube->textureCoord(0,1); cube->position(-0.5,-0.5,0.0);cube->normal(-0.408248,-0.816497,-0.408248);cube->textureCoord(1,1); cube->position(-0.5,-0.5,1.0);cube->normal(-0.666667,-0.333333,0.666667);cube->textureCoord(1,0); cube->position(0.5,-0.5,0.0);cube->normal(0.666667,-0.333333,-0.666667);cube->textureCoord(0,1); cube->position(0.5,0.5,0.0);cube->normal(0.408248,0.816497,-0.408248);cube->textureCoord(1,1); cube->position(0.5,-0.5,1.0);cube->normal(0.408248,-0.816497,0.408248);cube->textureCoord(0,0); cube->position(0.5,-0.5,0.0);cube->normal(0.666667,-0.333333,-0.666667);cube->textureCoord(1,0); cube->position(-0.5,-0.5,0.0);cube->normal(-0.408248,-0.816497,-0.408248);cube->textureCoord(0,0); cube->position(-0.5,0.5,1.0);cube->normal(-0.408248,0.816497,0.408248);cube->textureCoord(1,0); cube->position(0.5,0.5,0.0);cube->normal(0.408248,0.816497,-0.408248);cube->textureCoord(0,1); cube->position(-0.5,0.5,0.0);cube->normal(-0.666667,0.333333,-0.666667);cube->textureCoord(1,1); cube->position(0.5,0.5,1.0);cube->normal(0.666667,0.333333,0.666667);cube->textureCoord(0,0); cube->triangle(0,1,2);      cube->triangle(3,1,0); cube->triangle(4,5,6);      cube->triangle(4,7,5); cube->triangle(8,9,10);      cube->triangle(10,7,8); cube->triangle(4,11,12);   cube->triangle(4,13,11); cube->triangle(14,8,12);   cube->triangle(14,15,8); cube->triangle(16,17,18);   cube->triangle(16,19,17); cube->end(); return cube;}ManualObject* PhysicsEntity::createPlane(Ogre::String name, Ogre::String matName){ ManualObject* plane = new ManualObject(name); plane->begin(matName); plane->position(-0.5,0,-0.5);plane->normal(0.408248,-0.816497,0.408248);plane->textureCoord(1,0); plane->position(-0.5,0.0,0.5);plane->normal(-0.408248,-0.816497,-0.408248);plane->textureCoord(0,1); plane->position(0.5,0.0,0.5);plane->normal(0.666667,-0.333333,-0.666667);plane->textureCoord(1,1); plane->position(0.5,0.0,-0.5);plane->normal(-0.666667,-0.333333,0.666667);plane->textureCoord(0,0); plane->triangle(0,1,2);      plane->triangle(0,2,3); plane->end(); return plane;}

Last, but not least, is the main.cpp file, where everything gets blended:

main.cpp:

And this is it.Now...Let's get busy!

CODE NAME: Brachistochrona

2010-03-31T14:29:00.000-07:00

Today I'm going to reveal a classified archive. Or so it seems...

Physicians have been keeping this knowledge from us simple mortals for centuries...hehehe...
While trying to understand the essentials on variational mechanics, I have come into this article: http://jazz.openfun.org/wiki/Brachistochrone.
The problem of the brachistochrone dates back to Newton (17th century) and appears mentioned frequently as one of the problems that can be solved be means of this refined technique.

In fact, this problem appears to be the actual trigger for the modern physics (or better stated, its solution), and reveals a completely counterintuitive phenomenon:

The quickest path is not the straight line

Amazingly, the solution for a bead on a wire going in the least time from point A to point B:

Is the cycloid represented in the picture.

As standard human, one always tends to think that the shortest path would be the quickest one. It is wrong.
Analogously, when it comes to minimize the action (see previous posts on the subject), it is very hard to imagine that the stationary points would not come at flattening the curve for the action.
This is what makes so hard to understand the variational principles that rule over Lagrangian, Eulerian and Hamiltonian mechanics.

If one only had been told...

I have made a spreadsheet with openoffice (www.openoffice.org) where this can be numerically perceived also.

The spreadsheet is here.

I have made it with the help of a paper by Jozef Hanc: The original Euler's calculus-of-variations method: Key to Lagrangian mechanics for beginners.
I found it extremely appropriate for my case, and also absolutely clarifying. Thanks mr Hanc

Ilustrating Lagrangian Variables

2010-03-28T17:27:00.000-07:00

In the last post we had a tentative on Variational Mechanics which resulted very stimulating.

In order settle ideas down, I have left the subject cooling in the fridge and I have fiddled around with the Lego Mindstorms(c) robot we have at the RGEX research group, trying to make a positioning system that is accurate enough for our purposes.
Surely Jaume (my thesis director) will be very happy to hear about it.
I know this is not the exact place for it, but it's my blog and I do what I want to...hehehe...anyway, whoever wants to see what we do there, it is available at: www.rgex.blogspot.com.

And now, back to work:

The intention of this post it to illustrate some relationships that exist in Mechanics and that are repeated ad-vomitum in every reference one finds on Variational Mechanics.
The main problem I have encountered is that explanations mix numerical methods with physical concepts, and more often that seldom everything is written in the cryptic language of mathematicians. Of course, in complete absence of figures or anything of the like.

So, here we go...
With the help of some coding in C++ and QtOctave (https://forja.rediris.es/projects/csl-qtoctave/), I have made some graphics for the simulation of a moving particle.
I have used the same deductive principle I have seen in the papers, which is not intuitive at all. Hopefully some drawings will make it more clear.

A moving particle

Instead of trying to obtain the formula in a straightforward manner, the Least Action Principle suggest one to imagine an arbitrary particle that changes its position from one place to another. So, that's exactly what I have done.
The code iterates some 100 times and assigns random values to the new position of our particle. Really simple stuff.

The red dots represent all 100 different positions the particle goes when random increments in x and y are applied.
The initial point is at (10,10), and all the following are dQ further (dQ is divided in X and Y components, of course). dQ gets a random value for each iteration.
Once we have the new positions, it is straightforward to apply physics concepts to get the values of Velocity, Acceleration, Kinetic Energy and Potential Energy that have to be applied to get the particle there:

V=dQ/dt

A=dV/dt

Ep=m*A*dQ

Ec=1/2(m*|V|²)

With a dt value of 0.1 and a mass m of 5, it is also a piece of cake to obtain a table with the respective values.

Work and energy

The magnitudes Work and Energy are directly related to the Ep and Ec values.
From a physical point of view, there are some very interesting things that can be observed when putting their values against the |dQ| values:

Distance |dQ| vs kynetic energy Ec

Distance |dQ| vs potential energy Ep

From the same experiment, the random values of distance show a correlation of two types with the respective energies:

Kynetic energy is linearly related to the distance covered by the particle. Its values depend solely on the module of the velocity parameter, which itself is linearly dependent on the distance.
Potential energy shows a more quadratic behaviour, provided its acceleration component is multiplied by the covered distance.

It is also interesting to have a look at the graph Ec vs Ep, where also the quadratic behaviour is observed:

Kynetic energy Ec vs Potential energy Ep

Interestingly enough, both (|dQ| vs Ep) and (Ec vs Ep) are exactly the same shape, but for the scale of the X axis, which is of course linearly related.

The basic variational principle

The following step once these observations are made, is that of nesting our little freely moving particle loop into another broader loop, and accumulate the values of the position to form a trajectory.

Provided the inner loop generates a bunch of 100 probable future positions, one tempting criterion would be to choose that of minimum value for the total calculated energy...

The random accumulated trajectory of the particle

Just to verify that by adding the minimum energy substeps we still get a random trajectory.

Obviously, no constraints have been applied, and the particle represents a "Markov process".

What the basic variational principle teachs us is that, given certain constraints, the trajectory of the particle would have been that of minimum energy.

In our code for each time step the chosen subpath has been that of minimum energy.
Nonetheless, for the global trajectory, i.e. after a series of timesteps, a more holistic approach has to be taken in order to obtain the actual "natural" path.
That is the ODEs approach, and we'll talk about it in further episodes...

The following is the C++ code employed to obtain the tables. These have been plotted with QTOctave, an opensource version of MathLab.
Please note that it is the version which provides the complete path. The inner loop provides the random new positions and the outer loops concatenates one position after another, selecting each time that with the minimum energy.

main.cpp:

#include <stdio.h>#include <stdlib.h>#include <time.h>#include <math.h>using namespace std;int main(){ double Ec; double Ep; double Emin; double m; double dX; double dY; double dT; double Px=10; double Py=10; double Vx=0; double Vy=0; double VminX=0; double VminY=0; double dVx; double dVy; double Ax; double Ay; double Fx; double Fy; FILE * pFile; pFile = fopen ("myfile.txt","w"); dT=0.1; m=5; for (int i=0;i<1000;i++) {  for (int j=0;j<100;j++)  {   srand((unsigned) time(NULL));   dX=drand48() -drand48()  ;   dY=drand48() -drand48()  ;      printf("%f %f\n",dX,dY);   dVx=dX/dT-VminX;   dVy=dY/dT-VminY;   Vx=dX/dT;   Vy=dY/dT;   Ax=dVx/dT;   Ay=-9.8+dVy/dT;   Fx=m*Ax;   Fy=m*Ay;   Ec = 0.5*m*(sqrt(Vx*Vx + Vy*Vy));   Ep = Fx*dX + Fx*dX + Fy*dX + Fy*dY;   //   fprintf (pFile, "%f %f %f %f %f\n",Px,Py,Vx,Vy,Emin);   if (Ep+Ec<=Emin||j==0)   {    Emin=Ep+Ec;    Px=Px+dX;    Py=Py+dY;    VminX=Vx;    VminY=Vy;//    fprintf (pFile, "%f %f %f %f %f\n",Px,Py,Vx,Vy,Emin);   }  }  fprintf (pFile, "%f %f %f %f %f %f %f %f %f\n",Px,Py,sqrt(dX*dX+dY*dY),Vx,Vy,Ax,Ay,Ec,Ep); }}

ACTION!

2010-03-20T21:29:00.000-07:00

One more week spent on Lagrangian Dynamics...and more results!

Well, the main point about Lagrange and his ubiqutous formulae (that following level of abstraction I talked about in the last post), passes through the concept of ACTION.

Inexplicably, when learning physics at the Secondary or even at the Graduate studies, everything seems to finish when we get to the three laws of Newton.

That's the highest point one ever gets...and then we reach the concept of energy by integrating forces along paths.

But there is still one step further (or maybe more, now that the world seems to be round instead of flat, one never knows...).

The notion of action is that of the integral of the energy along time. Easy, isn't it? Its units are energy·time (Joule-second in SI).

Ok then. The problem, I guess, is that in order to be able to do something useful with it, one has to go also one step further in the mathematical notion of function: that of Functional.

And perhaps here is where our poor brains become at risk of melting...or maybe not, because a functional is just a function where the variable is...a function.

So, we are supposed to understand and memorize the linearity property of the integrals, and not the idea of functional? Come on...

Anyway, let's forgive our poor teachers of calculus and physics, who did their best I am sure, and get to work...

The discipline in charge to explain functionals is Variational Calculus, a fairly complex matter. When applied to physics, it happens to produce the Variational Principle of Physics, which is on where all Lagrangian, Jacobian, Hamiltonian and Eulerian dynamics rest.

So , in order to properly grasp the ideas of Lagrangian Dynamics, we better start with a little revision. The best one I found is from Mr. Claude Lacoursière's thesis "Ghost and machines: Regularized Variational Method for Interactive Simulation Multibodyes with dry friction" , combined with the following extremely clear explanation of the concept of action (http://wapedia.mobi/en/Action_(physics)):

1.-Work and Energy:

When a force is applied over a mass along a trajectory kinetic energy is dissipated. Integrating the formula, we can obtain its value as:

T(x')=(m·||x'||²)/2

where x' (aka velocity) is the derivative of the position, m is the mass of the particle and T(x') is the total kinetic energy.

For the cases when the force assumes values tangent to the energy function, then another form of energy appears: this is potential energy.

This special form of energy has the special property of being independent of the trajectory.

And this form of energy must also be considered into our energy equation if we want to account for everything.

This leaves our equation for the energy of the system the following way:

E=T(x')+V(x)

Being T(x') the kinetic energy described above and V(x) the potential energy.

2.-Basic Variational principle:

Once we have our Energy Equation, we are able to define the Action S for that Energy.
This would be, as explained above, the integral of this equation over time.

Here, t1 and t2 are the initial and final time positions, and L(x,x') is our LAGRANGIAN, which is exactly our Energy Equation, but with another name.

The funny thing about this integral is that it solves nothing by itself. There is a lot of possible solutions for our trajectory, all of them valid in mathematical terms, but only one of them is good for us.

The D'Alembert principle states that, "for a physical trajectory and infinitesimal displacements thereof at each point in time, and compatible with all imposed constraints, the work done by the virtual displacements along the trajectory vanishes when the motion is not bounded, and non-positive otherwise".

In plain English: the path followed is going to be that for which the action is minimized (or more strictly, for which the action is stationary, i.e. there is a minimum point).

To express the constraints mentioned by monsieur D'Alembert, all we need is to apply Euler-Lagrange equation:

∂L/∂x_i - d/dt ( ∂L/∂x_i' ) = 0

And then use the mathematical technique of the Lagrangian Multipliers (which I will not explain here), in order to simplify our Action Integral into something we can operate with, i.e. a second order differential equation.

And this is it. The basis, the cornerstone of Lagrangian Dynamics, demystified.
Once these concepts are comprehended, it is relatively easy to toddle around all the rest of this vast world...now at least we have legs!

LAGRANGIAN DYNAMICS

2010-03-14T19:02:00.000-07:00

These days I've been researching on Lagrangian Dynamics.

Finally, some light has been shed through the thick mesh of formulae, technical and mathematical verbiage, and an atom of understanding has entered into my banged head.

The paper that has finally put everything together is this:

http://box2d.googlecode.com/files/GDC2009_ErinCatto.zip

But before, I have had to go through an odissey of other papers to understand why they call it Lagrangian Dynamics, and Lagrangian by extension, to everything related to dynamics in the latest physics simulation R+D world.

The thing is as follows, explained in plane English for plane human beings (like myself, ahem...):

Classical Newtonian mechanics have some limitations when it comes to computing the equations of motion.
We want to iteratively solve a set of Ordinary Differential Equations (ODE) for acceleration, velocity and position if we want to know where our physical elements will be.
But it happened that mysteriously, when computing them, the total energy of the system got affected, because dissipative forces like damping, friction, etc... had to be also implemented, at that is a daunting task...

And there is where Lagrange comes into play: he theorized that we can make use of a further abstraction of our system and embed it into another, a bit broader conceptually, where all these forces would not be so determinant and we still could retrieve our results for position, by applying energetic conservation concepts.
Of course this abstraction comes at the price of having to understand a lot of many other things, not simple either (one funny thing about all this is that the only examples I have been able to find are that of a pendulum, a spring tied to a mass and a bead on a wire...If you didn't get it at the beginning, you'll never do...until you get it. Then these are so obvious you can't think of another!)

Anyway, the main point on Lagrangian approach to dynamics is that we can (and we normally do) find a set of equations that constraint the movement so the energy of the particles individually remains untouched.
These constraints materialize in the form of Jacobian Matrices or Lagrange Multipliers applied to each particle, rigid body or mass in a mass-spring system into our simulation (this is much deeper explained in the paper).

Fortunatelly also, all the above can be applied to Finite Element Method meshes. All we have to consider is that, in Lagrangian dynamical FEM, instead of having only one term for the stiffness matrix as usual:

[F]=[K]·[u]

a few more matrices that represent the dynamic terms appear:

[M][u''] + [C][u'] + [K][u] = [F]

And this is it regarding brand new concepts. It is as far as I can go by now.

In the way I have encountered many intensely related diverse topics, which should be explained in following posts. Some of them are:

Linear Complementary Problem Solvers
Principle of Virtual Work
Differential Equations
Smoothed Particle Hydrodynamics (SPH)
Kynematics

But for now I will leave this introduction as it is

Stiffness Method and SBRA

2010-03-10T19:51:00.000-08:00

Eventually I have managed to implement the full code for a simple Direct Stiffness calculator.
The previous post has been conveniently modified to reflect this.

The results have been fairly discouraging for my intemptions to use it as a basis for the SBRA: initial research led me to the wrong idea that the most computationally intensive part was that of assembly of the stiffness matrix, but my own experiments have proven the opposite.

The code provided, on a 1.66 GHz processor (the computer is actually dual core, but no parallelization was used), with 2GB RAM, takes 46 seconds only in the solving part. The modelled structure is 1000 nodes and 1800 elements large.

According to the paper, about 10⁷ iterations are needed over the solver to get a proper reliability assessment. Hence, grosso modo, applying SBRA to this structure on this computer would take 46*10⁷ seconds, some 14.5 long years.

It is likely that, for smaller structures and using parallelization, times become significantly smaller.

Nonetheless, the point I wanted to make here is that of the solving step being the most computationally intensive, instead of the presumed assembly step.

DSF Implementation

2010-03-05T04:50:00.000-08:00

Well...the following is the C++ code to achieve the initial stages of the DSF (e.g. member formation, globalization, merging and linear solving):
main.cpp:

#include <iostream>#include <math.h>#include "Stiffness.h"#include <vector>#include <gsl/gsl_matrix.h>#include <gsl/gsl_blas.h>#include <time.h>#include <ctime>#include <stdio.h>#include <unistd.h>#include <stdlib.h>int main(){ struct timespec start, stop1, stop2,stop3;    double accum;    clock_gettime( CLOCK_REALTIME, &start);    vector<Element> Elements;    vector<Node> Nodes;    Element Ele;    Node Node;    int numelements=0;    int numnodes=0;    //Create an structure of ten flat meshes one on top of the next    //In total makes 1000 nodes and 1800 elements    for (int h=0;h<10;h++)    {     for (int i=0;i<=9;i++)         {          for (int j=0;j<9;j++)          {           Ele.A=1;           Ele.E=100;           Ele.nodei=h*100+10*i+j;           Ele.nodej=h*100+10*i+j+1;           Elements.push_back(Ele);           numelements++;          }         }                  for (int i=0;i<=9;i++)         {          for (int j=0;j<9;j++)          {           Ele.A=1;           Ele.E=100;           Ele.nodei=h*100+10*j+i;           Ele.nodej=h*100+10*j+i+10;           Elements.push_back(Ele);           numelements++;          }         }                  for (int i=0;i<10;i++)            {              for (int j=0;j<10;j++)              {                  Node.Coordinates[0]=10*i;                  Node.Coordinates[1]=10*j;                  Node.Coordinates[2]=h*50;                  Node.DOF[0]=1;                  Node.DOF[1]=1;                  Node.DOF[2]=1;                  Node.Forces[0]=0;                  Node.Forces[1]=0;                  Node.Forces[2]=-10;                  Node.PrescribedDisplacements[0]=0;                  Node.PrescribedDisplacements[1]=0;                  Node.PrescribedDisplacements[2]=0;                  Nodes.push_back(Node);                                    numnodes++;              }            }    }        printf("Number of elements:%d\n",numelements);       printf("Number of nodes:%d\n",numnodes);     printf("Assembling...\n"); Structure *Fem=new Structure(Elements,Nodes); Fem->SpaceTrussMasterStiffness(); printf("Assembled!\n"); // printf("***********************************************************************\n");// printf("*******************************Kg Matrix*******************************\n");// printf("***********************************************************************\n");// for (int i=0;i<=Fem->NumNodes*3-1;i++)// {//  for (int j=0;j<=Fem->NumNodes*3-1;j++)//  {//   printf("%6.2f,",gsl_matrix_get(Fem->Kg,i,j));//  }//  printf("\n");// } clock_gettime( CLOCK_REALTIME, &stop1);    accum = ( stop1.tv_sec - start.tv_sec );    printf( "Elapsed time (ms): %4.6f\n", accum ); Fem->ModifiedMasterStiffness(Fem->Kg);//// printf("***********************************************************************\n");// printf("*************************Forces Vector*********************************\n");// printf("***********************************************************************\n");// for (int i=0;i<Fem->NumNodes*3;i++)// {//  printf("%g ", gsl_vector_get(Fem->VectorForces,i));// } Fem->ModifiedNodeForces();  printf("Solving...\n"); Fem->LinearSolve(Fem->Kg,Fem->VectorForces); delete Fem;    clock_gettime( CLOCK_REALTIME, &stop2);    accum = ( stop2.tv_sec - stop1.tv_sec );    printf( "Elapsed time (ms): %4.6f\n", accum );        clock_gettime( CLOCK_REALTIME, &stop3);    accum = ( stop3.tv_sec - start.tv_sec );    printf( "Total elapsed time (ms): %4.6f\n", accum );};

Element.h:
Element.cpp:
Node.h:
Node.cpp:

Stiffness.h:

Stiffness.cpp:

#include <iostream>#include <string>#include <cstdio>#include <vector>#include <math.h>#include <gsl/gsl_matrix.h>#include <gsl/gsl_blas.h>#include <gsl/gsl_linalg.h>#include "Stiffness.h"#include "Node.h"#include "Element.h"using namespace std;Structure::Structure(vector<Element> ImpElements,vector<Node> ImpNodes){ printf("Initializing...\n"); Elements=ImpElements; Nodes=ImpNodes; NumElements=Elements.size(); NumNodes=Nodes.size(); Kg = gsl_matrix_alloc(NumNodes*3,NumNodes*3); VectorForces=gsl_vector_alloc(NumNodes*3); int counter=0; int n1=Nodes.size(); int n2=3; for (int i=0;i<n1;i++) {  for (int j=0;j<n2;j++)  {   gsl_vector_set(VectorForces,counter,Nodes[i].Forces.at(j));   counter+=1;  } }// printf("Creating DOF table...\n"); this->PrescDisplacementDOFS();// printf("Created!\n");// printf("Initialized!\n");}Structure::~Structure(){ printf("Removing data...\n"); gsl_matrix_free( this->Kg );}void Structure::SpaceTrussMasterStiffness(){ int ii,jj; printf("Assembling...\n"); for (int e=0; e<NumElements; e++) {  gsl_matrix *Ke=gsl_matrix_alloc(6,6);  Element Ele=Elements.at(e);//  printf("**********************************************\n");//  printf("Processing element (%d)\n", e);//  printf("**********************************************\n");  int ni=Ele.nodei;  int nj=Ele.nodej;  int eftab[6]={3*(ni+1)-3,3*(ni+1)-2,3*(ni+1)-1,3*(nj+1)-3,3*(nj+1)-2,3*nj+1)-1};  SpaceBar2Stiffness(Nodes[ni].Coordinates,Nodes[nj].Coordinates,Ele.E,Ele.A,Ke);  int neldof=Ke->size1;  for (int i=0; i<=neldof-1; i++)  {   ii=eftab[i];   for (int j=i; j<=neldof-1; j++)   {    jj=eftab[j];    gsl_matrix_set(this->Kg,ii,jj,    gsl_matrix_get(this->Kg,ii,jj)    +gsl_matrix_get(Ke,i,j));    gsl_matrix_set(this->Kg,jj,ii,    gsl_matrix_get(this->Kg,ii,jj));   }  }  gsl_matrix_free( Ke );//  printf("Inserted!\n"); }}void Structure::ModifiedMasterStiffness(gsl_matrix *K){ int i,n,np; np=PrescDOFTable.size(); n=K->size1; printf("Modifying Kg...\n"); for (int k=0;k<np;k++) {  i=PrescDOFTable.at(k);  for (int j=0;j<n;j++)  {   gsl_matrix_set(K,i,j,0);   gsl_matrix_set(K,j,i,0);   gsl_matrix_set(K,i,i,1);  } }}void Structure::ModifiedNodeForces(){ int i,np,d,n,counter(0); printf("Modifying force vector...\n"); np=PrescDOFTable.size(); n=Kg->size1; int n1,n2; n1=Nodes.size(); n2=3; for (int i=0;i<n1;i++) {  for (int j=0;j<n2;j++)  {   gsl_vector_set(VectorForces,counter,Nodes[i].Forces.at(j));   counter+=1;  } } int c[n]; for (int k=0;k<np;k++) {  i=PrescDOFTable.at(k);  c[i]=0; } for (int k=0;k<np;k++) {  i=PrescDOFTable.at(k);  d=PrescForcesTable.at(k);  gsl_vector_set(VectorForces,i,d);  if (d==0)  {   continue;  }  for (int j=0;j<n;j++)  {   gsl_vector_set(VectorForces,j,gsl_matrix_get(Kg,i,j)*c[j]*d);  } } printf("\n");}void Structure::LinearSolve (gsl_matrix *a, gsl_vector *b){ gsl_vector *x=gsl_vector_alloc (a->size1); int s; gsl_permutation *p=gsl_permutation_alloc (a->size1); gsl_linalg_LU_decomp (a, p, &s); gsl_linalg_LU_solve (a, p, b, x); gsl_permutation_free (p);}void Structure::SpaceBar2Stiffness(vector<double> coor1,       vector<double> coor2,       double Em,       double A,       gsl_matrix *Ke){ double x21,y21,z21,EA,L,LL; x21=coor2.at(0)-coor1.at(0); y21=coor2.at(1)-coor1.at(1); z21=coor2.at(2)-coor1.at(2); EA=Em*A; LL=x21*x21+y21*y21+z21*z21; L=sqrt(LL); double Kaux[6][6] = {    { x21*x21, x21*y21, x21*z21,-x21*x21,-x21*y21,-x21*z21},    { y21*x21, y21*y21, y21*z21,-y21*x21,-y21*y21,-y21*z21},    { z21*x21, z21*y21, z21*z21,-z21*x21,-z21*y21,-z21*z21},    {-x21*x21,-x21*y21,-x21*z21, x21*x21, x21*y21, x21*z21},    {-y21*x21,-y21*y21,-y21*z21, y21*x21, y21*y21, y21*z21},    {-z21*x21,-z21*y21,-z21*z21, z21*x21, z21*y21, z21*z21}    }; for (int i=0; i<6;i++) {  for (int j=0; j<6;j++)   {    gsl_matrix_set(Ke,i,j,(EA/(L*L*L))* Kaux[i][j]);   } }}void Structure::PrescDisplacementDOFS(){ int k,index; k=0; for (int n=0;n<NumNodes;n++) {  for (int j=0;j<3;j++)  {   index=Nodes.at(n).DOF[j];   if (index>0)   {    PrescDOFTable.push_back(k+j);    PrescForcesTable.push_back(Nodes.at(n).Forces[j]);   }  }  k+=3; }}

I have used the GSL library for everything related to maths (http://www.gnu.org/software/gsl/).
This will make my life a lot easier once I manage to master it...Is fairly well documented, so it won't take long (I hope).
Please note that to make the timer work, under Ubuntu 9.10, I've had to activate the -lrt tag in the linker.

Direct Stiffness Method

2010-03-03T16:51:00.000-08:00

Well...as promised, these past days I've been researching on the possiblity of efficiently use the SBRA through the application of MonteCarlo methods.

I have been successful in my research and I have found this very interesting course:

http://www.colorado.edu/engineering/cas/courses.d/IFEM.d/

It is extremely complete!!

I have used it to refresh my rusted matrix knowledge from my university years, and so far I have managed to implement a very simple program - though effective for my purposes - that will allow me to benchmark the most time-consuming tasks in classical analysis.

This is a very brief resume of the algorithm from the course outlined above:

DIRECT STIFFNESS METHOD STEPS

Breakdown

Disconnection
Localization
Member Formation

Assembly

Globalization
Merge
Application of boundary conditions

Solution
Recovery of derived quantities

The code provided in the course is for Mathematica®, so all I have to do is translate it to C++.

I will provide the lines of code in the following posts.

Simulation Based Reliability Assesment (SBRA)

2010-02-25T11:42:00.000-08:00

This is the resume of the papers:

Both developed under supervision of Pavel Marek (see previous post about structural engineering conception).

The first one is very interesting since allows comparison between the actual "State of the Art" (i.e. Limit State approach to design), and the proposed explicitly probabilistic approach.
It describes a step-by-step process of analysis and design comparing both methods:

Loading and load combination

LRFD/Limit States

Each load is expressed by nominal value and its load factor
Load combinations are determined according to rules established in the Codes

SBRA

Each load is expressed by a load duration curve and its corresponding histogram
Load combination employs Monte Carlo sampling over some 10 million iterations of the analysis

Resistance and reference values

LRFD/Limit States

Resistances are combined with different resistance factors associated with failure mechanisms (yield stress, compression, flexure, shear,...) also provided by regulations
Design capacity=Nominal capacity x provided failure factor

SBRA

Reference value corresponds to the onset of the stress/deformation curve of the material when reaching yielding, or to a tolerable deformation
An histogram of yield stresses is used to feed each analysis iteration altogether with the histograms for the loads mentioned before

Frame analysis

LRFD/Limit States

One single iteration is made
Direct 2nd order analysis can be used, in order to check stabilities, but the most common approach is to do a 1st order analysis and then modificate it

SBRA

Analysis is repeated 10 million times using each time the values for loads and resistances extracted from the corresponding histograms. In these histograms, the more probable a load or resistance value, the more often will be used. An histogram is also called Probability Density Function in Monte Carlo literature.
Initial imperfections are explicitly taken into account
There is no need to check individual stability of columns (2nd order)
Resistance is related to the onset of yielding

Safety assessment

LRDF/Limit States

Once analysis is computed, for each element we check that the relationship Demand/Capacity is smaller than one

SBRA

Probability of failure Pf is compared against target probability Pd provided by the codes
P[(RV-LE)<0]=Pf

RV=yielding stress
LE=calculated stress

The proposal is very interesting in terms of providing a controlled design and an explicit probabilistic model of the structure.

Nevertheless, anyone would encounter the main drawback in the need for iterating 10 million times the whole structure, specially when it comes to real-time Lagrangian dynamics, as we are aimed to.
The steps for modelling the forces and displacements remain the same, it is, the traditional stiffness matrix.
Further research shows that the main time-consuming step into the stiffness matrix procedure (also applicable to FEM) is that of the assembly of the stiffness matrix (60-80%).
This means that, once a configuration is achieved and assembled, solving the matrix with different load cases should be relatively quick. Which is a relief if we want to implement SBRA.
Nevertheless, I still feel like I should pay deeper attention to this speed subject, and obtain my own benchmarking.

From Deterministic to Probabilistic Way of Thinking in Structural Engineering

2010-02-24T17:07:00.000-08:00

This is a very interesting article from the AECEF (Association of European Civil Engineering Faculties)!

From Deterministic to Probabilistic Way of Thinking in Structural Engineering
July 9, 2001. Adjusted Reprint: European Association of Civil Engineering Faculties - Newsletter 2001/ 1 p.

Prof. Ing. Pavel Marek, PhD., DrSc., F.ASCE, Prague pmarek@noise.cz
Prof. Jacques Brozzetti, Dr.h.c., Paris

Prestige of structural engineering
One of the speakers at the 1997 Congress of the American Society of Civil Engineers pointed out a significant drop of prestige of structural engineering. In the history of engineering education at the civil engineering departments in the USA this field used to occupy the highest echelon; therefore, the speaker was trying to find the explanation of this drop. He stated among others that the number of causes might include also the preference granted to other specializations due to general interests concerned with, for example, the development of transport networks and environment protection. In his opinion, however, the said drop might be connected also with the development of computers, structural design codes and corresponding software leading sometimes to such education of structural engineers which concentrates gradually on the application of sophisticated software, requiring from the designer merely the introduction of adequate input data, while the computer spouted almost immediately the dimensions of the structure, its reliability assessment and the whole documentation required by respective standards. In such process the designers need not even know the details of dimensioning or the substance of the reliability assessment. Is this opinion justified? Does the development of computer-aided structural design really belong among the causes of the above mentioned drop of prestige observed not only in the USA, but also in other countries?
Is the structural engineer the creator of structures or merely an interpreter of codes?
Together with the manufacturer and erector, the designer has been, and will always remain, the creator of the structure. His activities are based on professional knowledge, experience and cooperation with related professions. Although the computer revolution is providing ever more powerful instruments facilitating and accelerating his work, these instruments and how ever perfect codes can never replace the designer, responsible for effectiveness and reliability of his work.
The design codes and standards cannot cover all situations the designer may encounter, i.e., loading alternatives, performance conditions, etc. Often he has to decide himself on the basis of his own knowledge and experience in accordance with the „rules of the game“ of reliability assessment. In respect to safety, serviceability and durability of structures the development of codes and standards in the past few decades has resulted in a certain damping of the creative role of the designer who has become merely an interpreter of the rules and criteria formulated in the standard. The „rules of the game“ (i.e. the theoretical foundations of reliability assessment) of the Partial Factor Design (PFD, in the USA called Load and Resistance Factor Design, LRFD) are given in the codes excessively simplified and the designer during his education is not fully acquainted with their substance. The instructors often use the wording „the code states...“, thus avoiding the explanation of the problems with which they are often not thoroughly acquainted themselves. To be accurate, they cannot be thoroughly acquainted with them, as the commentaries and data explaining fully and consistently the background of the codes, various simplifications and the influence of calibration, are not available. The consequences of this development are that the students at departments of civil engineering are often educated primarily in the interpretation of codes and not in the creative engineering way of thinking. This fact must be afforded full attention.
Let us recall the related experience with the introduction of the PFD concept into AISC design codes in the USA In the field of steel structures a standard based on the LRFD method was issued in the USA. in 1986. This method is actually an analogy of the PFD method – see the Eurocodes. The code was introduced in order to replace the standards based on the deterministic Allowable Stress Design method (ASD). Although the issue of the LRFD standard was preceded by an extensive explanatory campaign and training courses emphasizing the advantages of the LRFD method, today – 15 years after it has been introduced – it is applied merely by one quarter of designers while the prevailing majority of designers in the country of the tallest buildings, biggest bridges and other unique structures, has remained faithful to the excessively simplified, but understandable deterministic ASD method (Iwankiw N. AISC, Chicago. Personal communication. 2000). This seemingly conservative attitude of the designers is usually explained by unsatisfactory teaching of the LRFD method at universities. However, the number of principal causes of reserve on the part of experienced USA designers may include their feeling that the LRFD method, developed in the pre-computer era, has been submitted to the designers too late and that it no longer provides qualitatively new possibilities corresponding with the computer era in respect of reliability assessment.
From slide-rule era to computer era in structural design
In the courses of steel, concrete and timber structures the students of civil engineering faculties may hear from some of their instructors that due to the introduction of the Eurocodes „nothing much will happen in the field of reliability assessment in the next few decades“. This is the statement which must be contradicted. The expansion of fast improving computers to the desk of every structural designer has produced profound qualitative and quantitative changes which have no analogy in the whole history of this field. The growing computer potential improves the prerequisites for the „re-engineering“ of the whole design process (i.e. its fundamental re-assessment and re-working) to adapt it to entirely new conditions and possibilities. We can follow with admiration the fast development of software for the analysis and dimensioning of structures according to existing codes and the production of their respective drawings. At the same time it is necessary to emphasize that entirely unsatisfactory attention has been afforded so far to the development of concepts and corresponding codes based on the qualitatively improved method of reliability assessment corresponding to computer potential available.
Since the early Sixties, many national and international codes for structural design based on deterministic concept have been replaced by a “semi-probabilistic” Partial Factors Design (PFD) such as that found in the Eurocodes. The PFD concepts have been developed using statistics, reliability theory and probability, however, without considering the computer revolution. The interpretation of the assessment format in codes is similar to the fully deterministic scheme applied in earlier codes except there are applied two partial factors instead of a single factor. The application of PFD does not require the designer to understand the rules hidden in the background of the codes. The semi-probabilistic background of the reliability assessment procedure has been considered by those writing the codes, however, the calibration and numerous simplifications introduced in the final format of codes affected the concept in such a way that the concept is better to be called “prescriptive” instead of “semi-probabilistic” (Iwankiw N. AISC, Chicago. Personal communication. 2000). The designer’s activities are limited to the interpretation of equations, criteria, instructions, factors, and „black boxes“ contained in the codes. The reliability check can be conducted using a calculator, slide rule or even long-hand, while the modern computer serves only as a “fast calculator”. The actual probability of failure and the reserves in bearing capacity cannot be explicitly evaluated using PFD codes. From the designer’s point of view, the application of PFD in practice is still deterministic. A designer‘s direct involvement in the assessment process is not assumed and, therefore, his/her creativity is suppressed.
Has the computer potential created the prerequisites for a qualitative improvement of the partial factors method? The answer can be illustrated by the following analogy: Is it sufficient to attach a high-efficient jet engine to the gondola of a balloon in order to achieve its incomparably higher velocity and efficiency? It can be concluded that the combination of the balloon and the jet engine will not create a higher quality means of transport. Analogously it is possible to conclude that the partial factors method based on numerous limitations and simplifying assumptions cannot be raised to a qualitatively higher level of the structural reliability assessment concept by its combination with computer potential. It can be concluded that the computer revolution leads to qualitatively new fully probabilistic structural reliability assessment concepts.
Application of elite research results to structural design codes
The results of elite research are usually applied to specific fields (offshore structures, space programs) by top-level experts. The conferences, however, lack the papers by research scientists explaining their ideas of the application of their concepts to codes and standards used by hundreds of thousands of designers in their everyday work. Who will bring the message from the elite researchers to the rank and file designers? An understandable explanation of scientific methods of reliability assessment used in the standards accepted by structural designers worldwide is a highly exacting task. However, without its solution the results of elite research remain merely the object of articles in scientific periodicals and the designer remains merely an interpreter of prescriptive codes.
Research affords attention to the development of risk engineering, fuzzy sets and other methods, while the designer lacks a fundamental, understandable and consistent method of determination of failure probability. Therefore, it is necessary to reassess the rules of the game of the reliability assessment, beginning with the load definition. The present day load expression in codes by the characteristic value and load factors must be replaced with a qualitatively higher form enabling to take into account also the loading history (such as, for example, the so called load duration curves, see book Simulation-based Reliability Assessment for Structural Engineers). With reference to bearing capacity it is necessary to provide a reference value applicable to the computation of the failure probability. Reliability should be expressed by a comparison of the computed failure probability with design target probability.
The preparedness of designer is the necessary prerequisite for the practical application of the probabilistic concept of reliability assessment. Let us turn our attention to the education of students at civil engineering departments and designers in post-graduate courses.
Deterministic or probabilistic approach in education?
Let us ask these questions: Is the approach applied in our courses to the solution of technical problems in structural design, deterministic or probabilistic? Are instructors infusing a deterministic understanding into the “knowledge-base” of their students, or is the fact that we are living in a world defined by random variables already accepted and applied in the educational process? In courses such as Statistics and Probability Models in Structural Engineering, the common textbooks are based on a “classical” approach to statistics and probability theory. Such an approach is limited to analytical and numerical solutions, and does not allow for transparent analysis of reliability functions that depend on the interaction of several random variables. The textbooks mostly remain silent on common real-world problems, such as the probability of failure of a structural component exposed to variable load combinations in which one might consider the contributions of variable yield stress, variable geometrical properties and random imperfections. In structural design courses, the interpretation and application of the existing codes are emphasized; however, students are using the codes without a full understanding of the actual reliability assessment rules and of the meanings of the factors used to express safety, durability, and serviceability of structural components.
Teaching reliability
Advances in computer technology allow for using simulation-based approach to solve numerous problems. The Monte Carlo simulation technique has been applied to basic problems in statistics and probability for a long time. The structural reliability assessment using direct Monte Carlo has been taught, for example, since 1989 at San Jose State University, California, and since 1996 at the Department of Civil Engineering, V©B TU Ostrava, Czech Republic, at the graduate and undergraduate levels. The positive response of the students, and their understanding encouraged the instructors. A team of undergraduate students developed, for example, a study proving that the LRFD method is not leading to a balanced safety (see Probabilistic Engineering Mechanics 14 (105 – 118), USA, 1998)). The new generation of civil engineers seems to be anxious to apply advanced computer technology to its fullest including application of simulation techniques in the analysis of multi-variable problems.
TeReCo project
With reference to the improvements expected in the field of structural reliability assessment the training of students and designers ranks among the most important tasks. What starting point should be chosen? A transition to the qualitatively higher probabilistic concepts will require the designer to change his way of thinking, i.e. to replace his current “deterministic thought-process” with the probabilistic one. The professional EC Committees consider the training of designers in this respect highly desirable. For this reason the Leonardo da Vinci Agency in Brussels has sponsored the TeReCo Project (Teaching Reliability Concepts using simulation techniques). This long term project had resulted in the work of 33 authors from nine countries being published in the textbook Probabilistic Assessment of Structures Using Monte Carlo Simulation. Background, Exercises and Software. (it is available since June 2001). The book acquaints the readers with the basis of a fully probabilistic structural reliability assessment concept, using as a tool the transparent SBRA method (Simulation-Based Reliability Assessment method documented in textbook Simulation-based Reliability Assessment for Structural Engineers and applied in about hundred papers). The concept allows for bypassing the „design-point” approach as well as the load and resistance factors, and leads to the reliability check expressed by P_f < P_d comparing the calculated probability of failure P_f with the target design probability P_d given in codes. The application of SBRA is explained in the book using 150 solved examples. On the attached CD-ROM, the reader will find the input files and computational tools enabling the duplication of the examples on a PC, the pilot data-base of mechanical properties (expressed by histograms) of selected structural steel grades, selected histograms (loads, imperfections, and more), manuals for computer programs and 55 selected presentations of examples (Microsoft PowerPoint). The book should serve as an aid in teaching undergraduate and graduate students and in introducing the designers to the strategy of the fully probabilistic reliability assessment of elements, components, members and simple systems using direct Monte Carlo simulation and modern PC computers.
Summary and conclusions
The structural engineering profession needs new approaches if we want to provide the best possible service to society. We have to consider the transition from the deterministic „way of thinking” to open-minded probabilistic concepts accepting the random character of individual variables involved as well as their interaction. Tools such as simulation techniques and powerful personal computers will contribute to reaching such goals. Students find these techniques easy to learn and thus they do not require the instructor to take a great deal of classroom time to explain the theoretical background. Once in the computer lab, students can explore to their hearts content and gain a fuller understanding of the effects of each parameter on the variability of the final answer. With this understanding students are better informed to make decisions about tradeoffs that need to be made, for example, between service life, durability and safety. The simulation technique should be included in the program of undergraduate and graduate studies and in corresponding textbooks to prepare students for the types of problems they will encounter in the real world, especially for the application of probabilistic structural reliability assessment concepts in the new generation of codes expected to be introduced in the near future. Such approach will make the engineer the creator of the structure and may bring the prestige of structural engineering back to one of the highest positions.

Limit State Design

2010-02-23T17:03:00.000-08:00

Limit State Design is the method now used in every regulation (Eurocode, Spanish CTE, USA LRFD, ...), so I think it is a primary target in the steepy road towards stochastic simulation.

Actually, a lot of bibliography points from Limit State Design to FORM (First Order Reliability Method).
Stochastic methods are being widely used into these environments (those of FORM).

But before, lets get a brief review on Limit States:

Limit States Design require the designer to establish a set of performance criteria (vibration levels, deflection, strength, stability, buckling, twisting, collapse,...) for the designed piece.
Then, according to its relevance, the variable of design enters within the area of the safety (ULS) or the comfort (SLS):

ULS (Ultimate Limit State): Satisfaction of ULS happens when all factored bending, shear, tensile and compressive stresses are below the factored resistances.

Stresses are magnified
Resistances are reduced

SLS (Serviceability Limit State): Satisfaction of SLS happens when deflection, vibration , fissuration,...are within a certain criteria

Hence, the actual method is then called, in European regulations, Partial Factors method, and in USA LRFD (Load and Resistance Factor Design).
Here the limits and the factors are prescripted by the regulatory organ, in a probabilistic manner, but are applied in the same deterministic way and within the same deterministic methodology as it has been for decades.

Collision Detection Summary

2010-02-20T21:08:00.000-08:00

Well, well...after a whole week of intense research, I think I already have some insight on Collision Detection so as to try to condense the matter in a post.

The following is a brief report made out of the references stated at the end:

INTRODUCTION

Collision detection is one of the major bottlenecks presented in any physics simulation given its computationally intensive idiosyncrasy.

It is also one of the main sources of error given that it serves as the basis for every other step in the simulation. Calculated collision points might later be treated as constraints within our solvers, or as force accumulators, and provide a lot of instabilities if they are not precise enough.

In general, collision detection / proximity queries / penetration depth algorithms are targeted to provide information about just that, to be used later for other things. These things are commonly physics simulations, haptic devices, or robotic simulations, and they require a real time speed of computation and robustness altogether.

A ROUGH REVIEW ON COLLISION DETECTION METHODOLOGIES FOR POLYGONAL MODELS

According to literature, it is possible to make an initial distinction depending on the nature of the graphic input data: polygonal models vs non-polygonal.

This separates the mostly common "soups of vertices" and triangles from solids, implicit surfaces and parametric surfaces.

Basically, all the researched collision detection methodologies focus on the first, as it is fairly easy to pre-process the others into polygon soups.

Furthermore, algorithms for non-polygonal surfaces are very specific and then would limit the input and the end user would be forced to provide information in this "format".

Once this distinction is made, each of the different algorithms focus on any or some of the following issues:

Pair processing vs N-processing: If a method emphasizes too much on improving precision in the measurement of the distance between two particular meshes, then it becomes too slow and so it has a limit to the amount of bodies the simulation can represent.
Static vs Dynamic: Others using the foreknown trajectories of the bodies require an amount of feedback between the dynamics engine and the collision detection package. Then the concept of frame coherence appears. These methods rely on the stable consistency of every time step to make calculations.
Rigid vs Deformable: When the meshes represent cloths or fluids (deformable meshes), as the dynamic characteristics of their shape is often complex and essentially non-linear, robustness and efficiency are frequently affected.
Large vs Small environment: There is also a problem generalizing the algorithms so they can serve as well in the short as in the long range, to avoid leaps in the numerical results that lead to instabilities and affect robustness.

BOUNDING VOLUME HIERARCHIES (BVHs)

An ubiquitous subject encountered in every researched method is that of the Bounding Volume Hierarchies.

These are a consequence of an iterative simplification of the problem of intersection of N bodies.

It consists in associating a simpler geometrical shape to our initial complex shape, so calculations of the distances can be made faster.

As long as couples of shapes gets closer and closer, we iteratively subdivide those rougher representations of them.
Many types have been devised:

k-DOP trees (Discrete Orientation Polytopes)
Octrees
R-trees
Cone trees
BSPs (Binary Space Partitions)

These trees are conformed of different simpler geometrical shapes that would be hierarchically sorted and that would represent our shapes within the tree.
There are various types:

Spheres
AABB (Axis Aligned Bounding Boxes)
OBBs (Oriented Bounding Boxes)
k-DOPs (Discrete Orientation Polytopes, where k means the number of faces)
SSV (Swept Sphere Volumes)

The main procedure consists of getting statistical data from all the vertices for each mesh in order to get the centroid (the mean value of all of them) and the maximum and minimum values to have the bounding volume (of whatever the chosen shape). Then it is possible to use the specific characteristics of the shape to proceed to the interference calculations.

With spheres, it is quite straightforward, having their radius and centre. With bounding boxes it is possible to "project" the mesh vertexes coordinates against their sides. With these projections, a set of lists is made for each dimension (called intervals) and then check interference sorting the lists.

Actually these trees, or hierarchies, are just an abstraction to get different levels of definition during computation.

When made between objects far apart, we speak about broadphase. In this state the nested resulting bounding volumes are simpler and rougher and have little to do with the actual shape they contain.

When objects get closer then appears the narrowphase. It triggers whenever an intersection between two bounding volumes has been detected in the broadphase. Implies the iterative subdivision of the bounding volume, either until a collision point is found, or a penetration distance is calculated, or a separation distance is obtained.

In this part many different approaches can be found, with very prolix performance results and assumptions.
Nevertheless, particularly one sounds louder than others: the GJK algorithm.

GJK ALGORITHM AND THE MINKOWSKY SUM

A principal character in this play is the GJK (Gilbert-Johnson-Keerthy) algorithm.

It implements an algebraic property of every pair of convex polyhedra: the Minkowsky sum.

Basically this sum (or difference, there is not much consensus but I am not a mathematician either), has the ability of returning a very important concept: when two convex polyhedra are added through it, if they share a part of the space (they intersect), the centroid of them is contained within the resulting polyhedron.

For every pair of objects that need to have their interference checked (remember, it triggers only between those objects that are close), all that has to be done is to iteratively see whether the features (vertices, edges and faces) are getting close or far from the centroid.
This is made using the concept of simplex, that can be applied to vertices, edges, triangles or tetrahedra (all of them known as polytopes).
Making a single loop that iterates through the mesh geometry, it is possible to discard the interference whenever a vector points out from the direction of the centroid of the Minkowsky sum. If none appears, then objects are overlapping.

This algorithm has an irrefutable efficiency, provided its popularity, and has proven to pass from an O(n2) complexity to O(nlogn) complexity.
Mainly it comes from the fact that gives a true / false answer very quickly, as it does not need to iterate the whole geometry for the answer (usually comes after a few iterations over the simplex checking algorithm).

The first problem one encounters is that distances are not straightforward in this method. The way to obtain penetration distance requires further processing and affects performance negatively when implemented. An error tolerance has to be provided and then successive iterations are made in a very similar manner to the "sweep and prune" way described above: testing that a "shrinking" polytope from the centroid is as small as this error.

Another drawback of this method is its limitation to convex polyhedra. Minkowsky sum property does not apply for concave polyhedra.
This can be tackled using a preprocessing convex shape decomposition that would provide us with sets of Bounding Volume Hierarchies nested within the initial concave shape.

Last but not least is the instability problem that arises when very small shapes are tested against larger ones. Due to the difference in size, the Minkowsky sum presents extremely oblong shaped facets. This commonly leads to infinite loops.

OPEN SOURCE COLLISION DETECTION

The following is a list describing all available Open Source algorithms I have tested and helped me out to figure out the described above:

I-COLLIDE collision detection system: is an interactive and exact collision-detection library for environments composed of many convex polyhedra or union of convex pieces, based on the expected constant time, incremental distance computation algorithm and algorithms to check for collision between multiple moving objects.
RAPID interference detection system: is a robust and accurate polygon interference detection library for pairs of unstructured polygonal models. It is applicable to polygon soups models which contain no adjacency information and obey no topological constraints. It is most suitable for close proximity configurations between highly tesselated smooth surfaces.
V-COLLIDE collision detection system: is a collision detection library for large dynamic environments, and unites the nbody processing algorithm of I-COLLIDE and the pair processing algorithm of RAPID. It is designed to operate on large numbers of static or moving polygonal objects to allow dynamic addition or deletion of objects between timesteps.
SOLID interference detection system: is a library for interference detection of multiple three-dimensional polygonal objects including polygon soups undergoing rigid motion. Its performance and applicability is comparable to that of V-COLLIDE.
SWIFT++: SWIFT++ is a collision detection package capable of detecting intersection, performing tolerance verification, computing approximate and exact distance, or determining the contacts between pairs of objects in a scene composed of general rigid polyhedral models. It is a major extension of the SWIFT system previously released from UNC. Uses the convex polyhedra decompostion.
OPCODE: It is similar to popular packages such as SOLID or RAPID, but more memory-friendly, and often faster. Stands for OPtimized COllision Detection.
PQP proximity query package: It pre-computes a hierarchical representation
of models using tightting oriented bounding box trees (OBBTrees). At runtime, the algorithm traverses two such trees and tests for overlaps between oriented bounding boxes based on a separating axis theorem, which takes less than 200 operations in practice.

Collision Detection Sources

2010-02-19T08:23:00.000-08:00

Well, the last two days I've been researching on the topic of Collision Detection (as planned).

After a short fight against OPCODE (http://www.codercorner.com/Opcode.htm), which stands for OPtimized COllision DEtection, I've been forced to abandon under mere technicalities with its implementation in C++.

Just couldn't make it compile under gcc and EasyEclipse...anyway, thanks to Pierre Terdiman (http://www.codercorner.com/blog/).

I have switched to PQP (http://gamma.cs.unc.edu/SSV/), which claims to be as fast and memory efficient as OPCODE, and I have already managed to make it work...now let's strip it off!

Working Programme

2010-02-18T07:12:00.000-08:00

Well, today is the beginning of the new approach towards the thesis!

According to my calendar, I still have 7,5 months here in Cercedilla (60km north of Madrid), in my parents' house, with nothing to worry about but "ora et labora"...

So, it seems reasonable to try to develop a full physics environment of my own, where at least I can control and understand everything that happens on it...maybe later on I can go back to ODE or Bullet!

According to documentation on Bullet and ODE, common physics engines implement the following:

Collision detection
Rigid body dynamics
Constraints
Soft body dynamics

The thing is that I also want to implement Finite Element Analysis and Stochastic environment...So my research timeline could be as follows:

February: Collision detection+Stochastic environment
March: Rigid body dynamics+Stochastic environment
April: Rigid body dynamics+Finite Element Analysis
May: Rigid body dynamics+Constraints+Finite Element Analysis
June: Rigid body dynamics+Soft body dynamics
July: Finite Element Analysis+Soft body dynamics
August: Finite Element Analysis
September: Stochastic environment

Of course, it is important to have in mind that everything doesn't have to be perfectly implemented. A rough and working sample of each item should be enough, as long as the whole system allows for further improvements.

I count heavily on the already developed code (Bullet and ODE are open source) and on the extensive documentation available (which sometimes is more a burden than a relief).

There are also two important deadlines:

June for the presentation of the thesis proposal to the UPC jury
October for the mobility grant to Ljubljana

Well...time will tell!

Bullet Physics Engine

2010-02-17T17:44:00.000-08:00

After a few months banging my head against Bullet (http://bulletphysics.org), which apparently could make my life a lot easier:

Bullet is a Collision Detection and Rigid Body Dynamics Library. The Library is Open Source and free for commercial use, under the ZLib license. This means you can use it in commercial games, even on next-generation consoles like Sony Playstation 3.

I have chosen to abandon this strategy for the time being given the following conclusions:

Physics engines in general try to be as general as possible, becoming extremely complex when it comes to learn how they work.
Physics engines are normally game-oriented, and although I haven't been able to prove it yet, their lack of precision might become an important issue when it comes to solve engineering problems.
There is no way around in deeply understanding the numerical methods they employ.

Nonetheless, this is the amount of research done so far:

Main available 3D physics engines are:

Open Source

Closed source/limited free distribution

Havok
Newton Game Dynamics
PhysX (formerly NovodeX and incorporating Meqon)
Vortex by CMLabs Simulations
Wild Pockets

The reason why I chose Bullet was mainly that it is open source. Secondly, when compared to the others (ODE and Tokamak), I found its basis and documentation a bit more solid. Thirdly, it was the only one I could put to work altogether with Ogre3D (http://www.ogre3d.org), the scene graph manager of my choice, where I can make simulations visible.
I have made a couple of videos available in youtube:

A set of bricks hit by a mass:

http://www.youtube.com/watch?v=tFJu2oo18Zs

A peristyle hit by a ball:

(Not yet available)

When trying to access the Softbody module of the engine, I have found it already too complex to link with Ogre3D.
Joints make the system unstable, with frequent segmentation faults and even automatic resets of my laptop.

With this in mind, I think the next step will be to focus on achieving the knowledge to tackle my own physics engine!

Thesis Proposal

2010-02-14T05:41:00.000-08:00

Introduction

State-of-the-art

Non-deterministic methods have been succesfully applied in science since the 50s to tackle uncertanty and to evaluate complexity.

Lately, these very methods, combined with modern calculus tools (Finite Elemet Method, Applied Element Method, Discrete Element Method, ...) are proving very helpful in automotive, aerospatial and naval industries to achieve sofisticated, reliable and precise designs.

On the other hand, new simulation tools have appeared in the last decade that were not possible before the increase in computer power and their clustering and parallelizing capabilities.

Videogame, film and virtual reality industries are using them and simulations flood the media with a quality that often makes it very difficult to separate fiction from truth.

In the building environment, perhaps given some inherent differences in the requirements of the product provided to society, we are still pretty far from those technological advances.

This could be seen as a disadvantegeous position but opens, in fact, a whole world of new opportunities for research.

The thesis proposed here aims to provide a frame where these two not so different disciplines (non-deterministic design methods and computer physics simulation) combine in order to open another path to better design tools.

The determinstic approach to structural design

Modern building design and analysis is almost exclusively based on stiffnes matrix method or, in best cases, on Finite Element software.

These softwares have reached today a considerable level of sophistication and versatility. However, one increasingly important aspect of analysis that these programs are unable to address is that of uncertainty in structural parameters and in loading and boundary conditions.

It is a well known fact that deterministic single- point evaluation of the response may under many circumstances produce an over-designed and excessively conservative system if the presence of parameter scatter is not taken into account.

International building codes, historically, have taken this procedure for granted and nowadays Limit States is the csompulsory method for evaluating any building (Eurocodes, ASCE, ACI,...).

These Limit States are given to the designers on a probability basis but have to be necessarily included into our deterministic analysis in the form of material strenght minoration and force majoration.

An stochastic approach allows us to define the reliability of our design in a different way: by calculating in an straighforward manner the probabilty of failure (or limit state achieving), and then comparing it against the standards.

Ordinary Differential Equations and physics simulation

The second weakness of this deterministic approach to the design of building structures is that tackling non-linearity (buckling of columns, terrain, earthquakes, …) has become an extremely unprecise field, plenty of approximate methods.

Intensive particle-based Lagrangian methods, on the other hand, is a relatively recent field of research (even though it uses classical newtonian physics), where the phenomena stated above simply arises as a consequence of the simultaneous behaviour.

Thesis targets

Achieve a computer tool with the following features:
Real-time ODE based physics.
Behavior-monitored structural elements and parameters.
Real-time design visualization and designer interaction.
Stochastic methods applied to different structural systems and probability-based evaluation of their reliability.
Building forensics of existing or failed buildings
Evaluate the adequateness of Lagrangian particle-based methods for engineering purposes.

http://www.mnbvlabs.com/Thesis/Thesis proposal.pdf