Element stiffness Matrix C + +

Source: Internet
Author: User
Tags modulus

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 + +

Related Article

Contact Us

The content source of this page is from Internet, which doesn't represent Alibaba Cloud's opinion; products and services mentioned on that page don't have any relationship with Alibaba Cloud. If the content of the page makes you feel confusing, please write us an email, we will handle the problem within 5 days after receiving your email.

If you find any instances of plagiarism from the community, please send an email to: info-contact@alibabacloud.com and provide relevant evidence. A staff member will contact you within 5 working days.

A Free Trial That Lets You Build Big!

Start building with 50+ products and up to 12 months usage for Elastic Compute Service

  • Sales Support

    1 on 1 presale consultation

  • After-Sales Support

    24/7 Technical Support 6 Free Tickets per Quarter Faster Response

  • Alibaba Cloud offers highly flexible support services tailored to meet your exact needs.