sábado, 26 de noviembre de 2011

Tetrahedra discretization of a beam

In this post I will introduce the advances of the last weeks regarding the discretization into tetrahedra of a beam with Grasshopper's VB component.

Rotations about an arbitrary axis
An important part of the research has been devoted to the topic of 3d rotation, as for locating the sub-nodes of the beam appropiately we need to transform the beam's local coordinates with respect of its axis into the actual world's coordinates for Rhino to interpret.
The most popular way to make it are rotation matrices (here there is a thorough explanation on how), but also quaternions have been contemplated and still not fully discarded. Apparently rotation matrices present some numerical drifting problems when utilized for many time steps, as well as other problems that quaternions solve quite easily (mentioned here).

Tetrahedra subdivision
The beam is defined by two points arbitrarily located in 3d space. These points define a director vector and also a distance.
The remaining geometric properties (height h and width b) in the attached code have been assimilated to a percentage of the lenght. Also the number of segments s has been hardcoded to be proportional to the width but it could be anything else.
The attached code uses two functions: one for coordinate transformation by means of transformation matrices and another for matrix multiplication. The rest of the subroutine simply locates points along the beam according to the width b, the height h and the segment separation s.

The image above shows how the beam is discretized into nodes located at the edges. Each segment gives place to five embedded tetrahedron (right).
Local axes of the beam will be later utilized for the matter integration.

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
    '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
      '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
      '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
      '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
    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))
      Line = New Rhino.Geometry.Line(nodes(i + 1), nodes(i + 2))
      Line = New Rhino.Geometry.Line(nodes(i + 2), nodes(i + 3))
      Line = New Rhino.Geometry.Line(nodes(i + 3), nodes(i))
      'Longitudinal segments of the edges
      Line = New Rhino.Geometry.Line(nodes(i), nodes(i + 4))
      Line = New Rhino.Geometry.Line(nodes(i + 1), nodes(i + 5))
      Line = New Rhino.Geometry.Line(nodes(i + 2), nodes(i + 6))
      Line = New Rhino.Geometry.Line(nodes(i + 3), nodes(i + 7))
    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))
      Line = New Rhino.Geometry.Line(nodes(i + 5), nodes(i + 8))
      'Right side faces
      Line = New Rhino.Geometry.Line(nodes(i + 2), nodes(i + 5))
      Line = New Rhino.Geometry.Line(nodes(i + 5), nodes(i + 10))
      'Bottom faces
      Line = New Rhino.Geometry.Line(nodes(i + 2), nodes(i + 3))
      Line = New Rhino.Geometry.Line(nodes(i + 3), nodes(i + 10))
      'Left side faces
      Line = New Rhino.Geometry.Line(nodes(i), nodes(i + 3))
      Line = New Rhino.Geometry.Line(nodes(i + 3), nodes(i + 8))

      'Sectional diagonals
      Line = New Rhino.Geometry.Line(nodes(i), nodes(i + 2))
      Line = New Rhino.Geometry.Line(nodes(i + 3), nodes(i + 5))
    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))
    Line = New Rhino.Geometry.Line(nodes(nodes.count - 3), nodes(nodes.count - 2))
    Line = New Rhino.Geometry.Line(nodes(nodes.count - 2), nodes(nodes.count - 1))
    Line = New Rhino.Geometry.Line(nodes(nodes.count - 1), nodes(nodes.count - 4))
    '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))
    '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
    'And then we return the resulting vector
    Return Vaux
  End Function

lunes, 31 de octubre de 2011

Yet one more article...

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

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

Thanks to my mentors, Jaume and Vojko

domingo, 9 de octubre de 2011

My first SCI Indexed Article!

Here it is! A whole achievement!

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

However, from all the article my favourite part is at the beginning: 
“Structural engineering is the art of molding materials we don’t wholly understand, into shapes we can’t fully analyze, so as to withstand forces we can’t really assess, in such a way that the community at large has no reason to suspect the extent of our ignorance.”
Which I found as a quote from James E. Amrhein of the Masonry Institute of America...

Well, happy autumn!

domingo, 25 de septiembre de 2011

First FEM + Lagrange multipliers program in Grasshopper

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

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

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

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

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

FEM Component:
LM Component:

LUPD Component:

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!

domingo, 17 de julio de 2011

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

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.


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.

lunes, 9 de mayo de 2011

Impulses, FEM and Verlet

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.

sábado, 19 de marzo de 2011


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

viernes, 4 de febrero de 2011


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.
    L = T - V.\quad
    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.