Stiffness Method and SBRA

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

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:
#ifndef ELEMENT_H_ #define ELEMENT_H_ class Element { public: Element(); virtual ~Element(); double A; double E; double L; int nodei; int nodej; double Sigma; double IntForces; }; #endif /*ELEMENT_H_*/ Element.cpp:
#include "Element.h" Element::Element() { } Element::~Element() { } Node.h:
#ifndef NODE_H_ #define NODE_H_ #include <iostream> #include <vector> using namespace std; class Node { public: Node(); ~Node(); vector<double> Coordinates; vector<double> DOF; vector<double> Forces; vector<double> PrescribedDisplacements; }; #endif /*NODE_H_*/ Node.cpp:
#include <iostream> #include "Node.h" using namespace std; Node::Node() { Coordinates.resize(3); DOF.resize(3); Forces.resize(3); PrescribedDisplacements.resize(3); } Node::~Node() { }  
Stiffness.h:
#ifndef STIFFNESS_H_ #define STIFFNESS_H_ #include <iostream> #include <vector> #include <gsl/gsl_matrix.h> #include <gsl/gsl_blas.h> #include "Node.h" #include "Element.h" using namespace std; class Structure { public: int NumNodes; int NumElements; vector<Node::Node> Nodes; vector<Element::Element> Elements; vector<double> PrescDOFTable; vector<double> PrescForcesTable; gsl_matrix *Kg; gsl_vector *VectorForces; Structure(vector<Element>,vector<Node>); ~Structure(); void SpaceTrussSolution(); void SpaceTrussMasterStiffness(); void ModifiedMasterStiffness(gsl_matrix*); void ModifiedNodeForces(); void SpaceTrussStresses(); void SpaceBar2Stiffness(vector<double>, vector<double>, double , double , gsl_matrix*); void PrescDisplacementDOFS(); void LinearSolve (gsl_matrix*,gsl_vector*); }; #endif /*Structure_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

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)

This is the resume of the papers:

• Safety assessment of a steel frame using LRFD and SBRA methods
• Probabilistic reliability assessment of a stell frame using the SBRA method

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.