## 位移法解平面桁架

1 1,3 i,j
2 1,2 i,j
3 3,4 i,j
4 2,4 i,j
5 3,2 i,j



 (1)

 (2)

 (3)

 (4)

 (5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

// ========================================================================= // Filename: Truss.cpp // Description: 位移法解平面桁架 // Version: 1.0 // Created: 2014-5-22 13:43:01 // Company: CCCC First Highway Consultants Co.,LTD // Revision: none // Compiler: g++ // Author: whchen (http://blog.whchen.net), // =========================================================================   #include #include #include #include #include #include #include #include #include&lt;Eigen/Dense&gt;   using namespace Eigen; using namespace std;   typedef struct node { double x,y; } node;   typedef struct frame { //iNum,jNum,E,A,L0,M 分别为节点号i,j,弹性模量,截面积,长度,杆单元刚度矩阵 int iNum,jNum; double E,A,L0; Matrix4d M; } frame;   typedef struct restrain { int num; int x,y; } restrain;   int main() { string sTmp; //IO ifstream fileInput("input.txt"); ofstream fileOutput("output.txt");   //变量定义 int NE,NJ,NP,NS;//单元数、节点数、节点荷载数,约束节点数   //输入单元节点荷载总数 getline(fileInput,sTmp); istringstream line(sTmp); line&gt;&gt;NE&gt;&gt;NJ&gt;&gt;NP&gt;&gt;NS;   node joint[NJ]; //输入节点信息 for(int i=0;i&lt;NJ;++i) { int k; getline(fileInput,sTmp); istringstream line(sTmp); line&gt;&gt;k; k-=1; line&gt;&gt;joint[k].x&gt;&gt;joint[k].y; }   frame elem[NE]; //输入单元信息 for(int i=0;i&lt;NE;++i) { int k,ii,jj; getline(fileInput,sTmp); istringstream line(sTmp); line&gt;&gt;k; k-=1; line&gt;&gt;ii&gt;&gt;jj&gt;&gt;elem[k].E&gt;&gt;elem[k].A; ii-=1; jj-=1; elem[k].iNum=ii; elem[k].jNum=jj; elem[k].L0=sqrt(pow(joint[jj].x-joint[ii].x,2)+pow(joint[jj].y-joint[ii].y,2));   double cosA,sinA; cosA=(joint[jj].x-joint[ii].x)/elem[k].L0; sinA=(joint[jj].y-joint[ii].y)/elem[k].L0;   //计算单元刚度矩阵 double aMat[4][4]={ cosA*cosA, sinA*cosA, -cosA*cosA, -sinA*cosA, cosA*sinA, sinA*sinA, -cosA*sinA, -sinA*sinA, -cosA*cosA, -sinA*cosA, cosA*cosA, sinA*cosA, -cosA*sinA, -sinA*sinA, cosA*sinA, sinA*sinA };   for(int i=0;i&lt;4;++i) for(int j=0;j&lt;4;++j) elem[k].M(i,j)=aMat[i][j];   elem[k].M*=(elem[k].E*elem[k].A/elem[k].L0); }   //叠加获得总刚度矩阵 MatrixXd KM=MatrixXd::Zero(2*NJ,2*NJ); for(int i=0;i&lt;NE;++i) { int ii=elem[i].iNum; int jj=elem[i].jNum;   KM(2*ii ,2*ii )+=elem[i].M(0,0); KM(2*ii ,2*ii+1)+=elem[i].M(0,1); KM(2*ii+1,2*ii )+=elem[i].M(1,0); KM(2*ii+1,2*ii+1)+=elem[i].M(1,1);   KM(2*ii ,2*jj )+=elem[i].M(0,2); KM(2*ii ,2*jj+1)+=elem[i].M(0,3); KM(2*ii+1,2*jj )+=elem[i].M(1,2); KM(2*ii+1,2*jj+1)+=elem[i].M(1,3);   KM(2*jj ,2*ii )+=elem[i].M(2,0); KM(2*jj ,2*ii+1)+=elem[i].M(2,1); KM(2*jj+1,2*ii )+=elem[i].M(3,0); KM(2*jj+1,2*ii+1)+=elem[i].M(3,1);   KM(2*jj ,2*jj )+=elem[i].M(2,2); KM(2*jj ,2*jj+1)+=elem[i].M(2,3); KM(2*jj+1,2*jj )+=elem[i].M(3,2); KM(2*jj+1,2*jj+1)+=elem[i].M(3,3); }   //定义总外力向量,位移向量 VectorXd VF=VectorXd::Zero(2*NJ); VectorXd VV=VectorXd::Zero(2*NJ); bool T[2*NJ]; memset(T,true,sizeof(T));   //输入荷载 for(int i=0;i&lt;NP;++i) { getline(fileInput,sTmp); istringstream line(sTmp); int num; double Px,Py; line&gt;&gt;num&gt;&gt;Px&gt;&gt;Py; num-=1; VF(2*num)=Px; VF(2*num+1)=Py; }   restrain rst[NS]; //输入约束信息 for(int i=0;i&lt;NS;++i) { getline(fileInput,sTmp); istringstream line(sTmp); line&gt;&gt;rst[i].num&gt;&gt;rst[i].x&gt;&gt;rst[i].y; rst[i].num-=1; VF(2*rst[i].num )=rst[i].x; if(rst[i].x) T[2*rst[i].num ]=false; VF(2*rst[i].num+1)=rst[i].y; if(rst[i].y) T[2*rst[i].num+1]=false; }   //开始写求解部分 //提取子矩阵，先求位移 int cnt=0; for(int i=0;i&lt;2*NJ;++i) if(T[i]) ++cnt; MatrixXd subM=MatrixXd::Zero(cnt,cnt); int ii,jj; ii=jj=0; for(int i=0;i&lt;2*NJ;++i) for(int j=0;j&lt;2*NJ;++j) if(T[i]&amp;&amp;T[j]) { subM(ii,jj)=KM(i,j); ++ii; if(ii&gt;=cnt) { ii=0; ++jj; } }   //提取子向量 VectorXd subVF=VectorXd::Zero(cnt); ii=0; for(int i=0;i&lt;2*NJ;++i) if(T[i]) subVF(ii++)=VF(i);   //求解位移 VectorXd subVV=subM.colPivHouseholderQr().solve(subVF);   //合并结果进总位移向量 ii=0; for(int i=0;i&lt;2*NJ;++i) if(T[i]) VV(i)=subVV(ii++);   VF=KM*VV;   //输出节点位移 fileOutput&lt;&lt;setw(25)&lt;&lt;"节点位移"&lt;&lt;endl; fileOutput&lt;&lt;" Joint"&lt;&lt;" X"&lt;&lt;" Y"&lt;&lt;endl; for(int i=0;i&lt;2*NJ;i+=2) { fileOutput&lt;&lt;setw(5)&lt;&lt;i/2+1&lt;&lt;" "; fileOutput&lt;&lt;fixed&lt;&lt;setprecision(12)&lt;&lt;setw(15)&lt;&lt;VV(i) &lt;&lt;" "&lt;&lt;setw(15)&lt;&lt;VV(i+1)&lt;&lt;endl; }   //输出单元轴力 fileOutput&lt;&lt;endl&lt;&lt;endl; fileOutput&lt;&lt;setw(25)&lt;&lt;"单元轴力"&lt;&lt;endl; fileOutput&lt;&lt;"Element"&lt;&lt;" F"&lt;&lt;endl; for(int i=0;i&lt;NE;++i) { int ii,jj; double E,A,L0; ii=elem[i].iNum; jj=elem[i].jNum; E=elem[i].E; A=elem[i].A; L0=elem[i].L0; double xi,xj,yi,yj,ui,uj,vi,vj; xi=joint[ii].x; yi=joint[ii].y; xj=joint[jj].x; yj=joint[jj].y; ui=VV(2*ii); vi=VV(2*ii+1); uj=VV(2*jj); vj=VV(2*jj+1); double fi=E*A/L0*((xj-xi)*(uj-ui)+(yj-yi)*(vj-vi))/L0; fileOutput&lt;&lt;setw(5)&lt;&lt;i+1&lt;&lt;" "; fileOutput&lt;&lt;fixed&lt;&lt;setprecision(12)&lt;&lt;setw(15)&lt;&lt;fi&lt;&lt;endl; }   //输出支座反力 fileOutput&lt;&lt;endl&lt;&lt;endl; fileOutput&lt;&lt;setw(25)&lt;&lt;"支承反力"&lt;&lt;endl; fileOutput&lt;&lt;" Joint"&lt;&lt;" FX"&lt;&lt;" FY"&lt;&lt;endl; for(int i=0;i&lt;NS;++i) { int num=rst[i].num; fileOutput&lt;&lt;setw(5)&lt;&lt;num+1&lt;&lt;" "; fileOutput&lt;&lt;fixed&lt;&lt;setprecision(12)&lt;&lt;setw(15)&lt;&lt;VF(2*num) &lt;&lt;" "&lt;&lt;setw(15)&lt;&lt;VF(2*num+1)&lt;&lt;endl; } return 0; }

单元数 节点数 集中荷载数 约束数量 节点1 x坐标 y坐标 节点2 x坐标 y坐标 节点3 x坐标 y坐标 节点4 x坐标 y坐标 节点5 x坐标 y坐标 ...... 单元1 i节点 j节点 弹性模量 截面积 单元2 i节点 j节点 弹性模量 截面积 单元3 i节点 j节点 弹性模量 截面积 单元4 i节点 j节点 弹性模量 截面积 单元5 i节点 j节点 弹性模量 截面积 ...... 荷载1 x分量 y分量 荷载2 x分量 y分量 荷载3 x分量 y分量 ...... 约束1 X是否固定 Y是否固定（1或0） 约束2 X是否固定 Y是否固定（1或0） 约束3 X是否固定 Y是否固定（1或0） ......

» 本博客采用署名 2.5 中国大陆许可协议进行许可，本文版权归作者所有，欢迎转载，但必须在明显位置给出原文连接。