矩阵位移法大作业监工日记

5/15

新建文件夹，先做好前期结构数据化的手写性梳理，写一点变量的定义，感觉有很多重复的信息，但是懒得精简了。

编写输入文件in.txt，和直接用const差不多，感觉不如直接程序内打表，因为后面结构数据化的时候还是直接存的信息，泛用化水平进一步降低。

  1 /* 矩阵位移法大作业
  2  单元数：10  结点数：12  整体结点位移编码：22 
  3  结点编码：
  4  1（0,0，1）
  5  2（2,3,4）
  6  3（5,6,7,）
  7  4（0,0,0）
  8  5（8,9,10）
  9  6（5,6,11）
 10  7（0,0,12）
 11  8（13,14,15）
 12  9（16,17,18）
 13  10（13,14,19）
 14  11（0,0,0） 
 15  12（20,21,22） 
 16  单元定位向量： 柱：0 梁：1
 17  （1） 【0,0,1,2,3,4】 1,2 0
 18  （2）【2,3,4,5,6,7】 2,3  1
 19  （3）【0,0,0,8,9,10】 4,5  0
 20  （4）【8,9,10,5,6,11】 5,6  0
 21  （5）【8,9,10,13,14,15】 5,8  1
 22  （6）【5，6,11,16,17,18】 6,9  1
 23  （7）【0,0,12,13,14,15】 7,8  0
 24  （8）【 13,14,15,16,17,18】 8,9  0
 25  （9）【13,14,19,20,21,22】 10,12  1
 26  （10）【0,0,0,20,21,22】 11,12  0
 27 */
 28 /* 输入文件 in.txt
 29  12 10 22
 30  100000 15000
 31  1000000 10000
 32  3 4 3
 33  4 3 2 
 34  20
 35  12 14 10
 36  0 0 1 
 37  2 3 4
 38  5 6 7
 39  0 0 0
 40  8 9 10
 41  5 6 11
 42  0 0 12
 43  13 14 15
 44  16 17 18
 45  13 14 19
 46  0 0 0
 47  20 21 22
 48  1 2 0
 49  2 3 1
 50  4 5 0
 51  5 6 0
 52  5 8 1
 53  6 9 1
 54  7 8 0
 55  8 9 0
 56  10 12 1
 57  11 12 0
 58 */ 
 59 
 60 
 61 
 62 
 63 #include<iostream>
 64 #include<fstream>
 65 #include<cstdio>
 66 #include<algorithm>
 67 #include<cmath> 
 68 using namespace std;
 69 double EA[3],EI[3]; //梁和柱的强度参数表
 70 double l[4],h[4],q[4],Fp; // 跨长、层高、均布荷载、集中力
 71 double K[25][25],k[25][25]; // 整体刚度矩阵 和 单元刚度矩阵
 72 double alpha,T[8][8]; //旋转角和旋转矩阵 
 73 int Njoint,Nelem,NglbDOF; //结点数，单元数，结点位移个数 
 74 struct Joint{ //结点信息 
 75     double x,y; //坐标
 76     int GDOF[5]; //结点位移编码 
 77 }J[20];
 78 struct Elem{ //单元信息 
 79     int N1,N2; //两个结点
 80     int GlbDOF[8]; //单元定位向量
 81      double len,A,EI,EA; //长度，转角，抗弯刚度，抗拉刚度 
 82 }E[15];
 83 int main(){
 84     cin>>Njoint>>Nelem>>NglbDOF; //读入信息 并 进行结构数据化 
 85     cin>>EA[0]>>EI[0]>>EA[1]>>EI[1];
 86     for(int i=1;i<=3;++i) cin>>l[i];
 87     for(int i=1;i<=3;++i) cin>>h[i]; cin>>Fp;
 88     for(int i=1;i<=3;++i) cin>>q[i];
 89     
 90     J[1].x=J[2].x=0;  //结点初始化 
 91     J[3].x=J[4].x=J[5].x=J[6].x=l[1];
 92     J[7].x=J[8].x=J[9].x=J[10].x=l[1]+l[2];
 93     J[11].x=J[12].x=l[1]+l[2]+l[3];
 94     J[1].y=J[4].y=J[7].y=J[11].y=0;
 95     J[12].y=h[3];
 96     J[5].y=J[8].y=J[10].y=h[2];
 97     J[2].y=h[1];
 98     J[3].y=J[6].y=J[9].y=2.0*h[2];    
 99     for(int i=1;i<=Njoint;++i) 
100         for(int j=1;j<=3;++j) 
101             cin>>J[i].GDOF[j];
102             
103     for(int i=1;i<=Nelem;++i){ //单元初始化 
104         cin>>E[i].N1>>E[i].N2;
105         int opt; cin>>opt;
106         E[i].EA=EA[opt]; E[i].EI=EI[opt];
107         for(int j=1;j<=3;++j){
108             E[i].GlbDOF[j]=J[E[i].N1].GDOF[j];
109             E[i].GlbDOF[j+3]=J[E[i].N2].GDOF[j];
110         } 
111         E[i].len=sqrt(pow(J[E[i].N1].x-J[E[i].N2].x,2)+pow(J[E[i].N1].y-J[E[i].N2].y,2));
112         E[i].A=
113     }
114 
115     
116 }

5/15

5/16

开始写单元刚度矩阵的生成，一个元素一个元素敲，接着生成旋转矩阵，开始进行整体刚度阵的集成。矩阵乘法直接抄的一年级时做的大作业，估计后面解线性方程组的时候还会再copy。

  1 /* 矩阵位移法大作业
  2  单元数：10  结点数：12  整体结点位移编码：22 
  3  结点编码：
  4  1（0,0，1）
  5  2（2,3,4）
  6  3（5,6,7,）
  7  4（0,0,0）
  8  5（8,9,10）
  9  6（5,6,11）
 10  7（0,0,12）
 11  8（13,14,15）
 12  9（16,17,18）
 13  10（13,14,19）
 14  11（0,0,0） 
 15  12（20,21,22） 
 16  单元定位向量： 柱：0 梁：1
 17  （1） 【0,0,1,2,3,4】 1,2 0
 18  （2）【2,3,4,5,6,7】 2,3  1
 19  （3）【0,0,0,8,9,10】 4,5  0
 20  （4）【8,9,10,5,6,11】 5,6  0
 21  （5）【8,9,10,13,14,15】 5,8  1
 22  （6）【5，6,11,16,17,18】 6,9  1
 23  （7）【0,0,12,13,14,15】 7,8  0
 24  （8）【 13,14,15,16,17,18】 8,9  0
 25  （9）【13,14,19,20,21,22】 10,12  1
 26  （10）【0,0,0,20,21,22】 11,12  0
 27 */
 28 /* 输入文件 in.txt
 29  12 10 22
 30  100000 15000
 31  1000000 10000
 32  3 4 3
 33  4 3 2 
 34  20
 35  12 14 10
 36  0 0 1 
 37  2 3 4
 38  5 6 7
 39  0 0 0
 40  8 9 10
 41  5 6 11
 42  0 0 12
 43  13 14 15
 44  16 17 18
 45  13 14 19
 46  0 0 0
 47  20 21 22
 48  1 2 0
 49  2 3 1
 50  4 5 0
 51  5 6 0
 52  5 8 1
 53  6 9 1
 54  7 8 0
 55  8 9 0
 56  10 12 1
 57  11 12 0
 58 */ 
 59 
 60 
 61 
 62 
 63 #include<iostream>
 64 #include<fstream>
 65 #include<cstdio>
 66 #include<algorithm>
 67 #include<cmath> 
 68 #include<string.h>
 69 using namespace std;
 70 double EA[3],EI[3]; //梁和柱的强度参数表
 71 double l[4],h[4],q[4],Fp; // 跨长、层高、均布荷载、集中力
 72 // 整体刚度矩阵 和 单元刚度矩阵
 73 //旋转角和旋转矩阵 
 74 int Njoint,Nelem,NglbDOF; //结点数，单元数，结点位移个数 
 75 struct Joint{ //结点信息 
 76     double x,y; //坐标
 77     int GDOF[5]; //结点位移编码 
 78 }J[20];
 79 struct Elem{ //单元信息 
 80     int N1,N2; //两个结点
 81     int GlbDOF[8]; //单元定位向量
 82      double len,A,EI,EA; //长度，转角，抗弯刚度，抗拉刚度 
 83 }E[15];
 84 struct Matrix{ //便于矩阵运算的方法 
 85     int n,m;
 86     double a[25][25];
 87     void clear(int _n=0,int _m=0){
 88         n=_n;m=_m; memset(a,0.0,sizeof(a));
 89     } 
 90     Matrix operator *(const Matrix& b)const {//矩阵乘法
 91         Matrix c;
 92         c.n=n; 
 93         c.m=b.m;
 94         int n=n, m=b.m, k=m;
 95         for (int i=1;i<=n;++i)
 96             for (int j=1;j<=m;++j)
 97                 for (int l=1;l<=k;++l)
 98                     c.a[i][j]+=a[i][l]*a[l][j];
 99         return c;
100     }
101 }K,_k,k,T,TT;// 整体刚度矩阵 和 单元刚度矩阵 旋转角和旋转矩阵
102 void init(){ //初始化过程 
103     cin>>Njoint>>Nelem>>NglbDOF; //读入信息 并 进行结构数据化 
104     cin>>EA[0]>>EI[0]>>EA[1]>>EI[1];
105     for(int i=1;i<=3;++i) cin>>l[i];
106     for(int i=1;i<=3;++i) cin>>h[i]; cin>>Fp;
107     for(int i=1;i<=3;++i) cin>>q[i];
108     
109     J[1].x=J[2].x=0;  //结点初始化 
110     J[3].x=J[4].x=J[5].x=J[6].x=l[1];
111     J[7].x=J[8].x=J[9].x=J[10].x=l[1]+l[2];
112     J[11].x=J[12].x=l[1]+l[2]+l[3];
113     J[1].y=J[4].y=J[7].y=J[11].y=0;
114     J[12].y=h[3];
115     J[5].y=J[8].y=J[10].y=h[2];
116     J[2].y=h[1];
117     J[3].y=J[6].y=J[9].y=2.0*h[2];    
118     for(int i=1;i<=Njoint;++i) 
119         for(int j=1;j<=3;++j) 
120             cin>>J[i].GDOF[j];
121             
122     for(int i=1;i<=Nelem;++i){ //单元初始化 
123         cin>>E[i].N1>>E[i].N2;
124         int opt; cin>>opt;
125         E[i].EA=EA[opt]; E[i].EI=EI[opt];
126         for(int j=1;j<=3;++j){
127             E[i].GlbDOF[j]=J[E[i].N1].GDOF[j];
128             E[i].GlbDOF[j+3]=J[E[i].N2].GDOF[j];
129         } 
130         E[i].len=sqrt(pow(J[E[i].N1].x-J[E[i].N2].x,2)+pow(J[E[i].N1].y-J[E[i].N2].y,2));
131         E[i].A=8.0*atan(1.0)-atan((J[E[i].N2].y-J[E[i].N1].y)/(J[E[i].N2].x-J[E[i].N1].x)); 
132     }
133     K.clear(22,22);_k.clear(6,6); k.clear(6,6);
134     T.clear(6,6);TT.clear(6,6);
135 }
136 void Integration(int i){
137     double _a=E[i].EA/E[i].len;
138     double _b=12.0*E[i].EI/pow(E[i].len,3);
139     double _c=6.0*E[i].EI/pow(E[i].len,2);
140     double _d=4.0*E[i].EI/E[i].len;
141     _k.a[1][1]=_k.a[4][4]=_a;
142     _k.a[1][4]=_k.a[4][1]=-_a;
143     _k.a[2][2]=_k.a[5][5]=_b;
144     _k.a[2][5]=_k.a[5][2]=-_b;
145     _k.a[2][3]=_k.a[3][2]=_k.a[2][6]=_k.a[6][2]=_c;
146     _k.a[3][5]=_k.a[5][3]=_k.a[6][5]=_k.a[5][6]=-_c;
147     _k.a[3][3]=_k.a[6][6]=_d;
148     _k.a[3][6]=_k.a[6][3]=0.5*_d; 
149     double _A=E[i].A; 
150     T.a[1][1]=T.a[2][2]=T.a[4][4]=T.a[5][5]=cos(_A);
151     T.a[1][2]=T.a[4][5]=sin(_A);
152     T.a[2][1]=T.a[5][4]=-sin(_A);
153     T.a[3][3]=T.a[6][6]=1.0;
154     for(int j=1;j<=6;++j)
155         for(int l=1;l<=6;++l)
156             TT.a[j][l]=T.a[l][j];
157     k=TT*_k*T;
158     for(int j=1;j<=6;++j)
159         for(int l=1;l<=6;++l){
160             int I=E[i].GlbDOF[j],J=E[i].GlbDOF[l];
161             if(I>0 && J>0)
162                 K.a[I][j]+=k.a[j][l];
163         }
164 }
165 void calcP(){
166     
167 }
168 int main(){
169     init(); //开始数据结构化
170     for(int i=1;i<=Nelem;++i)
171          Integration(i); //集成刚度阵  
172     calcP(); //计算整体结点荷载 
173 }

5/16

6/12

突然基本上写完了，完成了整体结点荷载向量的集成，从而计算出了位移和杆端力，完成了矩阵位移法的基本过程。调试过程主要是和各种除以0的精度问题做斗争，还发现大一时候写的高斯消元是错误的，有一处倒着枚举写成了正着枚举，而大作业自己给的例子过于简陋以至于没有发现这个问题，骗了郑老师的4.0

  1 /* 矩阵位移法大作业
  2  单元数：10  结点数：12  整体结点位移编码：22 
  3  结点编码：
  4  1（0,0，1）
  5  2（2,3,4）
  6  3（5,6,7,）
  7  4（0,0,0）
  8  5（8,9,10）
  9  6（5,6,11）
 10  7（0,0,12）
 11  8（13,14,15）
 12  9（16,17,18）
 13  10（13,14,19）
 14  11（0,0,0） 
 15  12（20,21,22） 
 16  单元定位向量： 柱：0 梁：1
 17  （1） 【0,0,1,2,3,4】 1,2 0
 18  （2）【2,3,4,5,6,7】 2,3  1
 19  （3）【0,0,0,8,9,10】 4,5  0
 20  （4）【8,9,10,5,6,11】 5,6  0
 21  （5）【8,9,10,13,14,15】 5,8  1
 22  （6）【5，6,11,16,17,18】 6,9  1
 23  （7）【0,0,12,13,14,15】 7,8  0
 24  （8）【 13,14,15,16,17,18】 8,9  0
 25  （9）【13,14,19,20,21,22】 10,12  1
 26  （10）【0,0,0,20,21,22】 11,12  0
 27 */
 28 /* 输入文件 in.txt
 29  12 10 22
 30  100000 15000
 31  1000000 10000
 32  3 4 3
 33  4 3 2 
 34  20
 35  12 14 10
 36  0 0 1 
 37  2 3 4
 38  5 6 7
 39  0 0 0
 40  8 9 10
 41  5 6 11
 42  0 0 12
 43  13 14 15
 44  16 17 18
 45  13 14 19
 46  0 0 0
 47  20 21 22
 48  1 2 0
 49  2 3 1
 50  4 5 0
 51  5 6 0
 52  5 8 1
 53  6 9 1
 54  7 8 0
 55  8 9 0
 56  10 12 1
 57  11 12 0
 58 */ 
 59 
 60 
 61 
 62 
 63 #include<iostream>
 64 #include<fstream>
 65 #include<cstdio>
 66 #include<algorithm>
 67 #include<cmath> 
 68 #include<string.h>
 69 using namespace std;
 70 const double eps=1e-18;
 71 double EA[3],EI[3]; //梁和柱的强度参数表
 72 double l[4],h[4],q[4],Fp; // 跨长、层高、均布荷载、集中力
 73 // 整体刚度矩阵 和 单元刚度矩阵
 74 //旋转角和旋转矩阵 
 75 int Njoint,Nelem,NglbDOF; //结点数，单元数，结点位移个数 
 76 struct Joint{ //结点信息 
 77     double x,y; //坐标
 78     int GDOF[5]; //结点位移编码 
 79 }J[20];
 80 struct Elem{ //单元信息 
 81     int N1,N2; //两个结点
 82     int GlbDOF[8]; //单元定位向量
 83      double len,A,EI,EA; //长度，转角，抗弯刚度，抗拉刚度 
 84 }E[15];
 85 struct Matrix{ //便于矩阵运算的方法 
 86     int n,m;
 87     double a[25][25];
 88     Matrix() { //默认构造函数
 89         n = 0; m = 0;
 90         memset(a, 0, sizeof(a));
 91     }
 92     void clear(int _n=0,int _m=0){
 93         n=_n;m=_m; memset(a,0.0,sizeof(a));
 94     } 
 95         Matrix operator +(const Matrix& b)const {//矩阵加法
 96         Matrix c;
 97         c.n=b.n; 
 98         c.m=b.m;
 99         int nn=c.n,mm=c.m;
100         for (int i=1;i<=nn;++i)
101             for (int j=1;j<=mm;++j)
102                 c.a[i][j]=a[i][j]+b.a[i][j];
103         return c;
104     }
105     Matrix operator *(const Matrix& b)const {//矩阵乘法
106         Matrix c;
107         c.n=n; 
108         c.m=b.m;
109         int nn=n,mm=b.m, kk=m;
110         for (int i=1;i<=nn;++i)
111             for (int j=1;j<=mm;++j)
112                 for (int l=1;l<=kk;++l)
113                     c.a[i][j]+=a[i][l]*b.a[l][j];
114         return c;
115     }
116     void Gauss() {//解线性方程组 
117         int r=0;
118     //    cout<<"n="<<n<<" m="<<m<<endl;
119         for (int i=1,j=1;i<=n&&j<=n;++i,++j){
120             int id=0;
121             for (int k=i;k<=n;++k)
122                 if (abs(a[k][j])>eps){
123                     id = k; break;
124                 }
125             if(!id) --i;
126             else{
127                 for(int k=j;k<=m;++k)
128                     swap(a[i][k],a[id][k]);
129             //    cout<<" look"<<endl;
130                 for(int k=m;k>=j;--k)
131                     a[i][k]/=a[i][j]/*,cout<<a[i][k]<<" "*/;
132                 //    cout<<endl;
133                 for (int k=i+1;k<=n;++k)
134                     if (abs(a[k][j])>eps){
135                         double K_=a[k][j]/a[i][j];
136                         for(int l=j;l<=m;++l){
137                             a[k][l]-= K_*a[i][l];
138                         }
139                     }
140                 //    else a[k][j]=0.0;
141                 ++r;
142             }
143     /*        cout<<"KKK"<<endl;
144         for(int i=1;i<=n;++i){
145             for(int j=1;j<=n+1;++j)
146             cout<<a[i][j]<<" ";
147             cout<<endl; 
148         }*/
149         }
150     /*    cout<<"KKK"<<endl;
151         for(int i=1;i<=n;++i){
152             for(int j=1;j<=n+1;++j)
153             cout<<a[i][j]<<" ";
154             cout<<endl; 
155         }*/
156         for(int i=n;i>=1;--i)
157             for(int j=i+1;j<=n;++j)
158                 a[i][n+1]-=a[j][n+1]*a[i][j];
159     }
160 /*    void out() const{//输出矩阵
161     cout<<"n="<<n<<" m="<<m<<endl;
162         for (int i = 1; i <= n; i++) {
163             for (int j = 1; j <= m; ++j)
164                 cout << a[i][j] << " ";
165             cout << endl;
166         }
167         cout << endl;
168     }*/
169 }K,_k,k[25],T,TT,P,PFe[25],Cal,Delta,delta[25];// 整体刚度矩阵  单元刚度矩阵 旋转角和旋转矩阵
170 void init(){ //初始化过程 
171     cin>>Njoint>>Nelem>>NglbDOF; //读入信息 并 进行结构数据化 
172     cin>>EA[0]>>EI[0]>>EA[1]>>EI[1];
173     for(int i=1;i<=3;++i) cin>>l[i];
174     for(int i=1;i<=3;++i) cin>>h[i]; cin>>Fp;
175     for(int i=1;i<=3;++i) cin>>q[i];
176     
177     J[1].x=J[2].x=0;  //结点初始化 
178     J[3].x=J[4].x=J[5].x=J[6].x=l[1];
179     J[7].x=J[8].x=J[9].x=J[10].x=l[1]+l[2];
180     J[11].x=J[12].x=l[1]+l[2]+l[3];
181     J[1].y=J[4].y=J[7].y=J[11].y=0;
182     J[12].y=h[3];
183     J[5].y=J[8].y=J[10].y=h[2];
184     J[2].y=h[1];
185     J[3].y=J[6].y=J[9].y=2.0*h[2];    
186     for(int i=1;i<=Njoint;++i) 
187         for(int j=1;j<=3;++j) 
188             cin>>J[i].GDOF[j];
189             
190     for(int i=1;i<=Nelem;++i){ //单元初始化 
191         cin>>E[i].N1>>E[i].N2;
192         int opt; cin>>opt;
193         E[i].EA=EA[opt]; E[i].EI=EI[opt];
194         for(int j=1;j<=3;++j){
195             E[i].GlbDOF[j]=J[E[i].N1].GDOF[j];
196             E[i].GlbDOF[j+3]=J[E[i].N2].GDOF[j];
197         } 
198         E[i].len=sqrt(pow(J[E[i].N1].x-J[E[i].N2].x,2)+pow(J[E[i].N1].y-J[E[i].N2].y,2));
199         double dy=J[E[i].N2].y-J[E[i].N1].y,dx=J[E[i].N2].x-J[E[i].N1].x;
200         if(abs(dx)<eps) dx=eps;
201         E[i].A=8.0*atan(1.0)-atan((dy)/(dx));
202     //    cout<<"atan(1)"<<8*atan(1.0)<<endl;
203     //    cout<<dy<<" "<<dx<<" "<<atan(dy/dx)<<endl;
204         k[i].clear(6,6);
205         PFe[i].clear(6,1);
206         delta[i].clear(6,1);
207     }
208     K.clear(NglbDOF,NglbDOF);_k.clear(6,6); 
209     T.clear(6,6);TT.clear(6,6);
210     P.clear(NglbDOF,1);
211     Cal.clear(NglbDOF,NglbDOF+1);
212     Delta.clear(NglbDOF,1);
213 }
214 void calT(double _A){
215     T.clear(6,6); TT.clear(6,6);
216 /*    for(int i=1;i<=6;++i){
217         for(int j=1;j<=6;++j){
218             cout<<TT.a[i][j]<<" ";
219         }cout<<endl;
220     }*/
221 //    cout<<"A="<<_A<<endl;
222     T.a[1][1]=T.a[2][2]=T.a[4][4]=T.a[5][5]=cos(_A);
223 //    cout<<"cos="<<cos(_A)<<endl;
224     T.a[1][2]=T.a[4][5]=sin(_A);
225     T.a[2][1]=T.a[5][4]=-sin(_A);
226     T.a[3][3]=T.a[6][6]=1.0;
227     for(int j=1;j<=6;++j)
228         for(int l=1;l<=6;++l)
229             TT.a[j][l]=T.a[l][j];
230 }
231 void Integration(int i){
232 //    cout<<"i="<<i<<endl;
233 //    cout<<"len="<<E[i].len<<endl;    
234     _k.clear(6,6);
235     double _a=E[i].EA/E[i].len;
236     double _b=12.0*E[i].EI/pow(E[i].len,3);
237     double _c=6.0*E[i].EI/pow(E[i].len,2);
238     double _d=4.0*E[i].EI/E[i].len;
239     _k.a[1][1]=_k.a[4][4]=_a;
240     _k.a[1][4]=_k.a[4][1]=-_a;
241     _k.a[2][2]=_k.a[5][5]=_b;
242     _k.a[2][5]=_k.a[5][2]=-_b;
243     _k.a[2][3]=_k.a[3][2]=_k.a[2][6]=_k.a[6][2]=_c;
244     _k.a[3][5]=_k.a[5][3]=_k.a[6][5]=_k.a[5][6]=-_c;
245     _k.a[3][3]=_k.a[6][6]=_d;
246     _k.a[3][6]=_k.a[6][3]=0.5*_d; 
247     double _A=E[i].A; 
248     calT(_A);
249 //    cout<<"fin"<<i<<endl;
250 /*    for(int i=1;i<=6;++i){
251         for(int j=1;j<=6;++j){
252             cout<<TT.a[i][j]<<" ";
253         }cout<<endl;
254     }
255         cout<<endl;
256     for(int i=1;i<=6;++i){
257         for(int j=1;j<=6;++j){
258             cout<<T.a[i][j]<<" ";
259         }cout<<endl;
260     }*/
261 //    TT.out();T.out();_k.out();
262 //    _k.out();
263 //    cout<<"chengfa"<<endl;
264 //    (TT*_k*T).out();
265     k[i]=(TT*_k*T);
266 //    k[i].out();
267     for(int j=1;j<=6;++j)
268         for(int l=1;l<=6;++l){
269             int I=E[i].GlbDOF[j],J=E[i].GlbDOF[l];
270             if(I>0 && J>0)
271                 K.a[I][J]+=k[i].a[j][l];
272         }
273 }
274 
275 
276 void calPE(int qi,int Ei){
277     double l=E[Ei].len,_A=E[Ei].A;
278 //    cout<<_A<<endl;
279     calT(_A);
280 //    cout<<"calt"<<endl;
281     PFe[Ei].clear(6,1);
282 //    cout<<"clear"<<endl;
283 //    cout<<"aaaa"<<endl;
284 //    cout<<"q="<<q[qi]<<endl;
285     PFe[Ei].a[1][1]=0.0;
286     PFe[Ei].a[2][1]=-q[qi]*l/2.0;
287     PFe[Ei].a[3][1]=-q[qi]*l*l/12.0; 
288     PFe[Ei].a[4][1]=0.0;
289     PFe[Ei].a[5][1]=-q[qi]*l/2.0;
290     PFe[Ei].a[6][1]=q[qi]*l*l/12.0; 
291 /*    for(int i=1;i<=6;++i){
292         for(int j=1;j<=6;++j){
293             cout<<TT.a[i][j]<<" ";
294         }cout<<endl;
295     }
296     cout<<"aa"<<TT.n<<" "<<TT.m<<endl;
297     cout<<"aa"<<PFe[Ei].n<<" "<<PFe[Ei].m<<endl; */
298 //    TT.out(); PFe[Ei].out();
299 //    (TT*PFe[Ei]).out();
300     PFe[Ei]=(TT*PFe[Ei]);
301 //    cout<<"PFE"<<endl;
302 //    PFe[Ei].out();
303     for(int i=1;i<=6;++i)
304         P.a[E[Ei].GlbDOF[i]][1]-=PFe[Ei].a[i][1];
305 //    cout<<"P"<<endl;
306 //    P.out();
307 }
308 void calcP(){
309     P.a[2][1]+=Fp;
310 //    cout<<"calp"<<endl;
311     calPE(1,5);calPE(2,6);calPE(3,9);
312 }
313 void caldelta(){
314     for(int i=1;i<=NglbDOF;++i){
315         for(int j=1;j<=NglbDOF;++j){
316             Cal.a[i][j]=K.a[i][j];
317         //    if(abs(Cal.a[i][j])<eps) Cal.a[i][j]=0.0;
318         }
319         Cal.a[i][NglbDOF+1]=P.a[i][1];
320     //    if(abs(Cal.a[i][NglbDOF+1])<eps) Cal.a[i][NglbDOF+1]=0.0;
321     }
322 //    cout<<"cal"<<endl;
323 //    Cal.out();
324 //    K.out();
325 //    P.out();
326     Cal.Gauss();
327 //    Cal.out();
328     for(int i=1;i<=NglbDOF;++i)
329         Delta.a[i][1]=Cal.a[i][NglbDOF+1];
330 //        cout<<"Delta"<<endl;
331 //        Delta.out();
332 //K.out(); P.out();
333 //    Delta=K*P;
334 //    Delta.out();
335 } 
336 void calcF(){
337 //    cout<<"begin"<<endl;
338     for(int i=1;i<=Nelem;++i){
339         for(int j=1;j<=6;++j)
340             delta[i].a[j][1]=Delta.a[E[i].GlbDOF[j]][1];
341 //        delta[i].out();
342     }
343     for(int i=1;i<=Nelem;++i){
344         PFe[i]=PFe[i]+k[i]*delta[i];
345 //        cout<<"i="<<i<<endl;
346 //        PFe[i].out(); 
347     }
348 } 
349 int main(){
350     init(); //开始数据结构化
351     //cout<<"11111"<<endl;
352     for(int i=1;i<=Nelem;++i)
353          Integration(i); //集成刚度阵  
354 //    cout<<"2222"<<endl;
355 //    K.out();
356     calcP(); //计算整体结点荷载 
357 //    cout<<"3333"<<endl;
358     caldelta();//计算整体结点位移 
359 //    cout<<"4444"<<endl;
360     calcF();//计算单元杆端力向量 
361 //    cout<<"5555"<<endl;
362     for(int i=1;i<=Nelem;++i){
363         cout<<"第"<<i<<"根杆件的全局坐标内力值"<<endl;
364         for(int j=1;j<=6;++j)
365             cout<<PFe[i].a[j][1]<<" ";
366         cout<<endl;
367     } 
368 //    Delta.out(); 
369     return 0;
370 }

6/12

6/24

卡着ddl开始做一些处理，绘制内力图需要单元杆端力，并编写程序文档。

#include<iostream>
#include<fstream>
#include<cstdio>
#include<algorithm>
#include<cmath> 
#include<string.h>
using namespace std;
const double eps=1e-12; //计算精度 
double EA[3],EI[3]; //梁和柱的强度参数表
double l[4],h[4],q[4],Fp; // 跨长、层高、均布荷载、集中力
int Njoint,Nelem,NglbDOF; //结点数，单元数，结点位移个数 
struct Joint{ //结点信息 
    double x,y; //坐标
    int GDOF[5]; //结点位移编码 
}J[20];
struct Elem{ //单元信息 
    int N1,N2; //两个结点
    int GlbDOF[8]; //单元定位向量
    double len,A,EI,EA; //长度，转角，抗弯刚度，抗拉刚度 
}E[15];
struct Matrix{ //便于矩阵运算的方法 
    int n,m;
    double a[25][25];
    Matrix() { //默认构造函数
        n = 0; m = 0;
        memset(a, 0, sizeof(a));
    }
    void clear(int _n=0,int _m=0){ //初始化函数 
        n=_n;m=_m; memset(a,0.0,sizeof(a));
    } 
    Matrix operator +(const Matrix& b)const {//矩阵加法
        Matrix c;
        c.n=b.n; 
        c.m=b.m;
        int nn=c.n,mm=c.m;
        for (int i=1;i<=nn;++i)
            for (int j=1;j<=mm;++j)
                c.a[i][j]=a[i][j]+b.a[i][j];
        return c;
    }
    Matrix operator *(const Matrix& b)const {//矩阵乘法
        Matrix c;
        c.n=n; 
        c.m=b.m;
        int nn=n,mm=b.m, kk=m;
        for (int i=1;i<=nn;++i)
            for (int j=1;j<=mm;++j)
                for (int l=1;l<=kk;++l)
                    c.a[i][j]+=a[i][l]*b.a[l][j];
        return c;
    }
    void Gauss() {//高斯消元法 解线性方程组 
        int r=0;
        for (int i=1,j=1;i<=n&&j<=n;++i,++j){
            int id=0;
            for (int k=i;k<=n;++k)
                if (abs(a[k][j])>eps){
                    id = k; break;
                }
            if(!id) --i;
            else{
                for(int k=j;k<=m;++k)
                    swap(a[i][k],a[id][k]);
                for(int k=m;k>=j;--k)
                    a[i][k]/=a[i][j];
                for (int k=i+1;k<=n;++k)
                    if (abs(a[k][j])>eps){
                        double K_=a[k][j]/a[i][j];
                        for(int l=j;l<=m;++l){
                            a[k][l]-= K_*a[i][l];
                        }
                    }
                ++r;
            }
        }
        for(int i=n;i>=1;--i)
            for(int j=i+1;j<=n;++j)
                a[i][n+1]-=a[j][n+1]*a[i][j];
    }
}K,_k,k[25],T,TT,P,PFe[25],Fe[25],Cal,Delta,delta[25];// 整体刚度矩阵  单元刚度矩阵 旋转矩阵（和逆矩阵） 整体结点荷载 全局坐标的单元杆端力 局部坐标的单元杆端力 用来gauss消元的矩阵 整体结点位移 每个杆件的结点位移 
void init(){ //初始化过程 
    cin>>Njoint>>Nelem>>NglbDOF; //读入信息 并 进行结构数据化 
    cin>>EA[0]>>EI[0]>>EA[1]>>EI[1];
    for(int i=1;i<=3;++i) cin>>l[i];
    for(int i=1;i<=3;++i) cin>>h[i]; cin>>Fp;
    for(int i=1;i<=3;++i) cin>>q[i];
    
    J[1].x=J[2].x=0;  //结点初始化 
    J[3].x=J[4].x=J[5].x=J[6].x=l[1];
    J[7].x=J[8].x=J[9].x=J[10].x=l[1]+l[2];
    J[11].x=J[12].x=l[1]+l[2]+l[3];
    J[1].y=J[4].y=J[7].y=J[11].y=0;
    J[12].y=h[3];
    J[5].y=J[8].y=J[10].y=h[2];
    J[2].y=h[1];
    J[3].y=J[6].y=J[9].y=2.0*h[2];    
    for(int i=1;i<=Njoint;++i) 
        for(int j=1;j<=3;++j) 
            cin>>J[i].GDOF[j];
            
    for(int i=1;i<=Nelem;++i){ //单元初始化 
        cin>>E[i].N1>>E[i].N2;
        int opt; cin>>opt;
        E[i].EA=EA[opt]; E[i].EI=EI[opt];
        for(int j=1;j<=3;++j){
            E[i].GlbDOF[j]=J[E[i].N1].GDOF[j];
            E[i].GlbDOF[j+3]=J[E[i].N2].GDOF[j];
        } 
        E[i].len=sqrt(pow(J[E[i].N1].x-J[E[i].N2].x,2)+pow(J[E[i].N1].y-J[E[i].N2].y,2));
        double dy=J[E[i].N2].y-J[E[i].N1].y,dx=J[E[i].N2].x-J[E[i].N1].x;
        if(abs(dx)<eps) dx=eps;
        E[i].A=8.0*atan(1.0)-atan((dy)/(dx));
        k[i].clear(6,6);
        PFe[i].clear(6,1);
        delta[i].clear(6,1);
    }
    K.clear(NglbDOF,NglbDOF);_k.clear(6,6); //各类矩阵/向量的初始化 
    T.clear(6,6);TT.clear(6,6);
    P.clear(NglbDOF,1);
    Cal.clear(NglbDOF,NglbDOF+1);
    Delta.clear(NglbDOF,1);
}
void calT(double _A){ //计算旋转矩阵 
    T.clear(6,6); TT.clear(6,6);
    T.a[1][1]=T.a[2][2]=T.a[4][4]=T.a[5][5]=cos(_A);
    T.a[1][2]=T.a[4][5]=sin(_A);
    T.a[2][1]=T.a[5][4]=-sin(_A);
    T.a[3][3]=T.a[6][6]=1.0;
    for(int j=1;j<=6;++j)
        for(int l=1;l<=6;++l)
            TT.a[j][l]=T.a[l][j];
}
void Integration(int i){ //集成刚度矩阵 
    _k.clear(6,6);
    double _a=E[i].EA/E[i].len;
    double _b=12.0*E[i].EI/pow(E[i].len,3);
    double _c=6.0*E[i].EI/pow(E[i].len,2);
    double _d=4.0*E[i].EI/E[i].len;
    _k.a[1][1]=_k.a[4][4]=_a;
    _k.a[1][4]=_k.a[4][1]=-_a;
    _k.a[2][2]=_k.a[5][5]=_b;
    _k.a[2][5]=_k.a[5][2]=-_b;
    _k.a[2][3]=_k.a[3][2]=_k.a[2][6]=_k.a[6][2]=_c;
    _k.a[3][5]=_k.a[5][3]=_k.a[6][5]=_k.a[5][6]=-_c;
    _k.a[3][3]=_k.a[6][6]=_d;
    _k.a[3][6]=_k.a[6][3]=0.5*_d; 
    double _A=E[i].A; 
    calT(_A);
    k[i]=(TT*_k*T);
    for(int j=1;j<=6;++j)
        for(int l=1;l<=6;++l){
            int I=E[i].GlbDOF[j],J=E[i].GlbDOF[l];
            if(I>0 && J>0)
                K.a[I][J]+=k[i].a[j][l];
        }
}

void calPE(int qi,int Ei){ //计算单元等效荷载 
    double l=E[Ei].len,_A=E[Ei].A;
    calT(_A);
    PFe[Ei].clear(6,1);
    PFe[Ei].a[1][1]=0.0;
    PFe[Ei].a[2][1]=-q[qi]*l/2.0;
    PFe[Ei].a[3][1]=-q[qi]*l*l/12.0; 
    PFe[Ei].a[4][1]=0.0;
    PFe[Ei].a[5][1]=-q[qi]*l/2.0;
    PFe[Ei].a[6][1]=q[qi]*l*l/12.0; 
    PFe[Ei]=(TT*PFe[Ei]);
    for(int i=1;i<=6;++i)
        P.a[E[Ei].GlbDOF[i]][1]-=PFe[Ei].a[i][1];
}
void calcP(){ //计算整体结点荷载 
    P.a[2][1]+=Fp;
    calPE(1,5);calPE(2,6);calPE(3,9);
}
void caldelta(){ //计算整体结点位移 
    for(int i=1;i<=NglbDOF;++i){
        for(int j=1;j<=NglbDOF;++j){
            Cal.a[i][j]=K.a[i][j];
        }
        Cal.a[i][NglbDOF+1]=P.a[i][1];
    }
    Cal.Gauss();
    for(int i=1;i<=NglbDOF;++i)
        Delta.a[i][1]=Cal.a[i][NglbDOF+1];
} 
void calcF(){ //计算杆端力向量 
    for(int i=1;i<=Nelem;++i){
        for(int j=1;j<=6;++j)
            delta[i].a[j][1]=Delta.a[E[i].GlbDOF[j]][1];
    }
    for(int i=1;i<=Nelem;++i){
        PFe[i]=PFe[i]+k[i]*delta[i];
        double _A=E[i].A;
        calT(_A);
        Fe[i]=T*PFe[i];
    }
} 
int main(){
    init(); //开始数据结构化
    for(int i=1;i<=Nelem;++i)
         Integration(i); //集成刚度阵  
    calcP(); //计算整体结点荷载 
    caldelta();//计算整体结点位移 
    calcF();//计算杆端力向量 
    for(int i=1;i<=Njoint;++i){ 
        cout<<"第"<<i<<"个结点的全局坐标位移值："<<endl;
        for(int j=1;j<=3;++j){
            if(abs(Delta.a[J[i].GDOF[j]][1])<eps) Delta.a[J[i].GDOF[j]][1]=0.0;
            cout<<Delta.a[J[i].GDOF[j]][1]<<" ";
        }
        cout<<endl;
    } 
    cout<<endl;
    for(int i=1;i<=Nelem;++i){
        cout<<"第"<<i<<"根杆件的全局坐标内力值："<<endl;
        for(int j=1;j<=6;++j){
            if(abs(PFe[i].a[j][1])<eps) PFe[i].a[j][1]=0.0; 
            cout<<PFe[i].a[j][1]<<" ";
        }
        cout<<endl;
    }
    cout<<endl; 
    for(int i=1;i<=Nelem;++i){
        cout<<"第"<<i<<"根杆件的单元杆端内力值："<<endl;
        for(int j=1;j<=6;++j){
            if(abs(Fe[i].a[j][1])<eps) Fe[i].a[j][1]=0.0; 
            cout<<(j<=2?-1:1)*Fe[i].a[j][1]<<" ";
        }
        cout<<endl;
    } 
    return 0;
} 
/* 矩阵位移法大作业
 单元数：10  结点数：12  整体结点位移编码：22 
 结点编码：
 1（0,0，1）
 2（2,3,4）
 3（5,6,7,）
 4（0,0,0）
 5（8,9,10）
 6（5,6,11）
 7（0,0,12）
 8（13,14,15）
 9（16,17,18）
 10（13,14,19）
 11（0,0,0） 
 12（20,21,22） 
 单元定位向量： 柱：0 梁：1
 （1） 【0,0,1,2,3,4】 1,2 0
 （2）【2,3,4,5,6,7】 2,3  1
 （3）【0,0,0,8,9,10】 4,5  0
 （4）【8,9,10,5,6,11】 5,6  0
 （5）【8,9,10,13,14,15】 5,8  1
 （6）【5，6,11,16,17,18】 6,9  1
 （7）【0,0,12,13,14,15】 7,8  0
 （8）【 13,14,15,16,17,18】 8,9  0
 （9）【13,14,19,20,21,22】 10,12  1
 （10）【0,0,0,20,21,22】 11,12  0
*/
/* 输入文件 in.txt
 12 10 22
 100000 15000
 1000000 10000
 3 4 3
 4 3 2 
 20
 12 14 10
 0 0 1 
 2 3 4
 5 6 7
 0 0 0
 8 9 10
 5 6 11
 0 0 12
 13 14 15
 16 17 18
 13 14 19
 0 0 0
 20 21 22
 1 2 0
 2 3 1
 4 5 0
 5 6 0
 5 8 1
 6 9 1
 7 8 0
 8 9 0
 10 12 1
 11 12 0
*/

6/24

posted @ 2024-05-15 23:26 Yu-shi 阅读(12) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

Yu-shi

西山苍苍

矩阵位移法大作业监工日记

公告