DSF Implementation

05 March 2010

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.

Rabindranat Subscriber|User role