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.