Element stiffness Matrix C + +

Because C + + does not encapsulate the Matrix class, it is still used in the computer commonly used numerical algorithm and program (c + +) a book of the header file "Matrix.h". The example in the book of finite element method is programmed, and the programming level is too much. The program is not so jumbled ...

Title: A three-strand four element triangle is given, and the modulus E and Poisson's ratio v are known and the unit thickness is t=1. The stiffness matrix of the unit is obtained;

Idea: According to the steps in the book, basically is the stiffness matrix ke=b ' *d*b*t*a, wherein, B is the element strain matrix, Db=s is the element stress matrix. To solve the problem step by step, the program is as follows:

#include <iostream>#include<cmath>#include"Matrix.h"//#include "Comm.h"using namespacestd;voidMainvoid){    DoubleAIJM (DoubleX2,DoubleX3,DoubleY2,Doubley3); DoubleBIJM (DoubleY2,Doubley3); DoubleCIJM (DoubleX2,DoubleX3);//Solving Parameters Ai,bj,am ... The function declaration    DoubleXi,xj,xm,yi,yj,ym;//i,j,m Point counter-clockwise arrangement    Doubleai,aj,am,bi,bj,bm,ci,cj,cm; //reading the three vertex coordinates of a triangular elementcout<<"xi,yi="; CIN>>xi>>Yi; cout<<"xj,yj="; CIN>>xj>>YJ; cout<<"xm,ym="; CIN>>xm>>ym; //solving parameter ai,aj,am,...Ai=AIJM (XJ,XM,YJ,YM); AJ=AIJM (XI,XM,YI,YM); Am=AIJM (XI,XJ,YI,YJ); Bi=BIJM (YJ,YM); BJ=BIJM (Ym,yi); BM=BIJM (YI,YJ); CI=CIJM (XJ,XM); CJ=CIJM (XI,XM); CM=CIJM (XI,XJ); //Output Parameterscout<<"ai="<<ai<<"\taj="<<aj<<"\tam="<<am<<Endl; cout<<"bi="<<bi<<"\tbj="<<bj<<"\tbm="<<bm<<Endl; cout<<"ci="<<ci<<"\tcj="<<cj<<"\tcm="<<cm<<Endl; /***********************************************/    DoubleE=2*pow (Ten,5);//modulus MPa    Doublev=0.2;//Poisson's ratio    Doublet=1;//unit thickness set to 1    DoubleA=3*4/2;//area of the triangular element    /***********************************************/    /*finding the strain matrix of a triangular element B=[BI,BJ,BM]*/    Const Doubleb[3][6]={{bi,0, BJ,0Bm0},        {0Ci0Cj0, CM}, {CI,BI,CJ,BJ,CM,BM}}; Matrix<Double> B (&b[0][0],3,6);//there is no const in front, otherwise tanspose transpose cannot be performedb=b/(2*A);//get strain matrix B;cout<<endl<<"b="<<Endl;                    Matrixlineprint (B); //output strain matrix B;    /*stress Matrix S=DB=[SI,SJ,SM];*/    Doublee0,v0,constant; //Plane Stress probleme0=E; V0=v; CONSTANT=e0/(2*a* (1-v0*v0));//Parametric constant    /************************************************* Const Double ssi[3][2]= {{Bi,v0*ci}, {V0*bi,ci},    {(1-v0) *CI/2, (1-v0) *BI/2}};    Const matrix<double> Si (&ssi[0][0],3,2);    cout<< "si=" <<endl;    Matrixlineprint (Si);    Const Double ssj[3][2]= {{BJ,V0*CJ}, {V0*BJ,CJ}, {(1-v0) *CJ/2, (1-v0) *BJ/2}};    Const matrix<double> SJ (&ssj[0][0],3,2);    cout<< "sj=" <<endl;    Matrixlineprint (SJ);    Const Double ssm[3][2]= {{bm,v0*cm}, {v0*bm,cm}, {(1-v0) *CM/2, (1-v0) *BM/2}};    Const matrix<double> SM (&ssm[0][0],3,2);    cout<< "sm=" <<endl;    Matrixlineprint (SM); ****************************************/    //synthetic stress matrix S;    Const Doubles[3][6]={{Bi,v0*ci,bj,v0*cj,bm,v0*cm}, {V0*bi,ci,v0*bj,cj,v0*bm,cm}, {(1-V0) *ci/2,(1-V0) *bi/2,(1-V0) *cj/2,(1-V0) *bj/2,(1-V0) *cm/2,(1-V0) *bm/2}    }; Matrix<Double> S (&s[0][0],3,6);//define the stress matrix S;S=s*constant;//get the stress matrix S;cout<<"s="<<Endl;                    Matrixlineprint (S); //Output Stress Matrix    /************** The transpose bt******************* of B strain Matrix*/Matrix<Double> BT (6,3);                Matrixtranspose (B,BT); //transpose The strain Matrix B to get BTcout<<"bt="<<Endl;                Matrixlineprint (BT); //output transpose strain matrix    /************** to find the element stiffness matrix ke***********************/Matrix<Double> Ke (6,6);//defining the element stiffness matrixKe= (bt*s) *t*a;//get the stiffness matrix of the unitcout<<"ke="<<Endl;                Matrixlineprint (Ke); //output single-rigid matrix}/*********************************************************/DoubleAIJM (DoubleX2,DoubleX3,DoubleY2,Doubley3) {    returnx2*y3-x3*y2;}/*********************************************************/DoubleBIJM (DoubleY2,Doubley3) {    returny2-Y3;}/*********************************************************/DoubleCIJM (DoubleX2,Doublex3) {    return-x2+x3;}

Input:

The output results are as follows:

Consistent with the results in the book, it is possible to verify the three properties of the element stiffness matrix: symmetry, singularity, main element greater than 0;

Element stiffness Matrix C + +