单元刚度矩阵 C++

　　由于C++没有封装矩阵类，所以还是用到了《计算机常用数值算法与程序》（C++）一书中的头文件“Matrix.h”。对《有限单元法》书中的例题进行了编程验算，编程水平太菜。程序冗杂得不行了...

　　题目：给出了一个勾三股四的单元三角形，弹性模量E和泊松比v已知，单元厚度t=1。求单元的刚度矩阵；

　　思路：按照书本中的步骤，基本是单刚矩阵Ke=B'*D*B*t*A;其中，B为单元应变矩阵，DB=S为单元应力矩阵。一步一步进行求解，程序如下：

#include<iostream>
#include<cmath>
#include"Matrix.h"
//#include"Comm.h"
using namespace std;
void main(void)
{
    double aijm(double x2,double x3,double y2,double y3);
    double bijm(double y2,double y3);
    double cijm(double x2,double x3);    //求解参数ai,bj,am...的函数声明

    double xi,xj,xm,yi,yj,ym;            //i，j，m点逆时针排列
    double ai,aj,am,bi,bj,bm,ci,cj,cm;
    //读取三角单元的三个顶点坐标
    cout<<"xi,yi=";
    cin>>xi>>yi;
    cout<<"xj,yj=";
    cin>>xj>>yj;
    cout<<"xm,ym=";
    cin>>xm>>ym;
    //求解参数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);
    //输出参数
    cout<<"ai="<<ai<<"\taj="<<aj<<"\tam="<<am<<endl;
    cout<<"bi="<<bi<<"\tbj="<<bj<<"\tbm="<<bm<<endl;
    cout<<"ci="<<ci<<"\tcj="<<cj<<"\tcm="<<cm<<endl;
    /***********************************************/
    double E=2*pow(10,5);   //弹性模量MPa
    double v=0.2;            //泊松比
    double t=1;                //单元厚度设为1
    double A=3*4/2;            //三角形单元的面积
    /***********************************************/
    /*求三角形单元的应变矩阵B=[Bi,Bj,Bm]*/
    const double b[3][6]=
    {
        {bi,0,bj,0,bm,0},
        {0,ci,0,cj,0,cm},
        {ci,bi,cj,bj,cm,bm}
    };
    matrix<double> B(&b[0][0],3,6);        //前面没有加const，否则后面不能进行tanspose转置
    B=B/(2*A);                            //得到应变矩阵B；
    cout<<endl<<"B="<<endl;
    MatrixLinePrint(B);                    //输出应变矩阵B；
    /*求应力矩阵S=DB=[si,sj,sm];*/
    double E0,v0,CONSTANT;
    //平面应力问题
    E0=E;
    v0=v;
    CONSTANT=E0/(2*A*(1-v0*v0));        //参量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);
    ****************************************/
    //合成应力矩阵S；
    const double s[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);        //定义应力矩阵S；
    S=S*CONSTANT;                        //得到应力矩阵S；
    cout<<"S="<<endl;
    MatrixLinePrint(S);                    //输出应力矩阵
    /**************求B应变矩阵的转置BT********************/
    matrix<double> BT(6,3);
    MatrixTranspose(B,BT);                //将应变矩阵B转置，得到BT
    cout<<"BT="<<endl;
    MatrixLinePrint(BT);                //输出转置应变矩阵
    /**************求单元刚度矩阵Ke***********************/
    matrix<double> Ke(6,6);                //定义单元刚度矩阵
    Ke=(BT*S)*t*A;                        //得到单元的刚度矩阵
    cout<<"Ke="<<endl;
    MatrixLinePrint(Ke);                //输出单刚矩阵


}
/*********************************************************/
double aijm(double x2,double x3,double y2,double y3)
{
    return x2*y3-x3*y2;
}
/*********************************************************/
double bijm(double y2,double y3)
{
    return y2-y3;
}
/*********************************************************/
double cijm(double x2,double x3)
{
    return -x2+x3;
}

输入：

输出结果如下：

对比书中的结果一致，可以验证单元刚度矩阵的三个性质：对称性，奇异性，主元素大于0;