-
Notifications
You must be signed in to change notification settings - Fork 1
Expand file tree
/
Copy pathmain.cpp
More file actions
200 lines (186 loc) · 7.48 KB
/
main.cpp
File metadata and controls
200 lines (186 loc) · 7.48 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
#include <cstdlib>
#include <string>
#include <iostream>
#include <math.h>
#include <fstream>
#include <iomanip>
#include "matrixops.h"
#include "node.h"
using namespace std;
int main()
{
ofstream Data;
Data.open ("Data.txt");
/*------------------Initialize Matrices, Nodes, and Variables-----------------*/
double**M = CreateMatrix(6); //Mass matrix
double**C = CreateMatrix(6); //Dampener matrix
double**K = CreateMatrix(6); //Spring matrix
double**summat = CreateMatrix(6); //Placeholder matrix for summing
double**summat2 = CreateMatrix(3); //Placeholder matrix for summing
double F[6]; //Force vector
double G[6]; //Gravitational force on all nodes
double a[6]; //Acceleration vector
double v[6]; //Velocity vector
double uc[6]; //Current displacement vector
double up[6]; //Previous displacement vector
double un[6]; //Next displacement vector
double sumvec[6]; //Placeholder vector for summing
double sumvec2[3]; //Placeholder vector for summing
double sumvec3[3];
double deltat = 0.0001; //Delta t
double deltat2 = deltat*deltat; //Delta t squared
Node n0; //Node 0
Node n1; //Node 1
Node n2; //Node 2
Node n3; //Node 3
Node n4; //Node 4
Node n5; //Node 5
/*------------------Set Matrices, Call initProb3 on Nodes---------------------*/
n0.initProb3(0);
n1.initProb3(1);
n2.initProb3(2);
n3.initProb3(3);
n4.initProb3(4);
n5.initProb3(5);
M[0][0] = n0.getselffactor(2); //Set mass matrix
M[1][1] = n1.getselffactor(2);
M[2][2] = n2.getselffactor(2);
M[3][3] = n3.getselffactor(2);
M[4][4] = n4.getselffactor(2);
M[5][5] = n5.getselffactor(2);
C[1][1] += n0.getconnfactor(0,1);
C[1][0] -= n0.getconnfactor(0,1);
C[0][1] -= n0.getconnfactor(0,1);
C[1][1] += n1.getconnfactor(0,1);
C[2][2] += n1.getconnfactor(0,1);
C[2][1] -= n1.getconnfactor(0,1);
C[1][2] -= n1.getconnfactor(0,1);
C[1][1] += n1.getconnfactor(1,1);
C[3][3] += n1.getconnfactor(1,1);
C[3][1] -= n1.getconnfactor(1,1);
C[1][3] -= n1.getconnfactor(1,1);
C[2][2] += n2.getconnfactor(0,1);
C[4][4] += n2.getconnfactor(0,1);
C[4][2] -= n2.getconnfactor(0,1);
C[2][4] -= n2.getconnfactor(0,1);
C[3][3] += n3.getconnfactor(0,1);
C[4][4] += n3.getconnfactor(0,1);
C[4][3] -= n3.getconnfactor(0,1);
C[3][4] -= n3.getconnfactor(0,1);
C[4][4] += n4.getconnfactor(0,1);
C[5][5] += n4.getconnfactor(0,1);
C[5][4] -= n4.getconnfactor(0,1);
C[4][5] -= n4.getconnfactor(0,1);
K[0][0] += n0.getconnfactor(0,0); //Set spring matrix
K[1][1] += n0.getconnfactor(0,0);
K[1][0] -= n0.getconnfactor(0,0);
K[0][1] -= n0.getconnfactor(0,0);
K[1][1] += n1.getconnfactor(0,0);
K[2][2] += n1.getconnfactor(0,0);
K[2][1] -= n1.getconnfactor(0,0);
K[1][2] -= n1.getconnfactor(0,0);
K[1][1] += n1.getconnfactor(1,0);
K[3][3] += n1.getconnfactor(1,0);
K[3][1] -= n1.getconnfactor(1,0);
K[1][3] -= n1.getconnfactor(1,0);
K[2][2] += n2.getconnfactor(0,0);
K[4][4] += n2.getconnfactor(0,0);
K[4][2] -= n2.getconnfactor(0,0);
K[2][4] -= n2.getconnfactor(0,0);
K[3][3] += n3.getconnfactor(0,0);
K[4][4] += n3.getconnfactor(0,0);
K[4][3] -= n3.getconnfactor(0,0);
K[3][4] -= n3.getconnfactor(0,0);
K[4][4] += n4.getconnfactor(0,0);
K[5][5] += n4.getconnfactor(0,0);
K[5][4] -= n4.getconnfactor(0,0);
K[4][5] -= n4.getconnfactor(0,0);
G[0] = n0.getselffactor(2)*(-9.80665);
G[1] = n1.getselffactor(2)*(-9.80665);
G[2] = n2.getselffactor(2)*(-9.80665);
G[3] = n3.getselffactor(2)*(-9.80665);
G[4] = n4.getselffactor(2)*(-9.80665);
G[5] = n5.getselffactor(2)*(-9.80665);
cout << "Mass matrix:" << endl;
MPrint(M,6,6);
cout << endl << "Dampener matrix:" << endl;
MPrint(C,6,6);
cout << endl << "Spring matrix:" << endl;
MPrint(K,6,6);
cout << endl;
for (int i=0; i<6; i++)
{
uc[i] = 0;
}
for (int i=0; i<6; i++)
{
up[i] = 0;
}
for (int i=0; i<6; i++)
{
un[i] = 0;
}
/*------------------Solution Start--------------------------------------------*/
Data << setw(8) << "0"; //Initial print to file
for (int printcount=0; printcount<6; printcount++)
{
Data << setw(18) << uc[printcount];
}
Data << endl;
for (double perm=0.0001; perm<=7.5; perm+=0.0001) //Permutations begin (0.0001s to 7.5s)
{
for (int i=0; i<6; i++)
{
F[i] = 0;
}
for (int i=0; i<6; i++)
{
F[i] += G[i];
}
subm2(K,M,summat,6,6,1,2/(deltat2));
mulsquvec(summat,uc,sumvec,6);
for (int i=0; i<6; i++)
{
F[i] -= sumvec[i];
}
subm2(M,C,summat,6,6,1/deltat2,(1/(2*deltat)));
mulsquvec(summat,up,sumvec,6);
for (int i=0; i<6; i++)
{
F[i] -= sumvec[i];
}
addm2(M,C,summat,6,6,1/deltat2,(1/(2*deltat)));
sumvec2[0] = F[1];
sumvec2[1] = F[2];
sumvec2[2] = F[4];
summat2[0][0] = summat[1][1];
summat2[0][1] = summat[1][2];
summat2[0][2] = summat[1][4];
summat2[1][0] = summat[2][1];
summat2[1][1] = summat[2][2];
summat2[1][2] = summat[2][4];
summat2[2][0] = summat[4][1];
summat2[2][1] = summat[4][2];
summat2[2][2] = summat[4][4];
lud(summat2,sumvec2,3,sumvec3);
un[1] = sumvec3[0];
un[2] = sumvec3[1];
un[4] = sumvec3[2];
for (int i=0; i<6; i++)
{
up[i] = uc[i];
}
for (int i=0; i<6; i++)
{
uc[i] = un[i];
}
cout << endl << perm;
Data << setw(8) << perm; //Print to file displacements
for (int printcount=0; printcount<6; printcount++)
{
Data << setw(18) << uc[printcount];
}
Data << endl;
}
system("pause");
}