-
Notifications
You must be signed in to change notification settings - Fork 0
/
main.cpp
139 lines (117 loc) · 3.17 KB
/
main.cpp
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
#include <iostream>
#include "Tools.h"
#include <math.h>
#include <Eigen/Eigen>
#define pi acos(-1)
using namespace std;
using namespace Eigen;
int main()
{
// physical domain
const double omega_l = 0.0;
const double omega_r = 2.0*pi;
const double g0 = Func_exact(omega_l); //Dirichlet BC
const double gL = Func_exact_x(omega_r); //Natural BC
// number of elements
const int nElem = 3;
// quadrature rule
const int nqp = 5;
double * qp = new double[nqp];
double * wq = new double[nqp];
Gauss(nqp, -1.0, 1.0, qp, wq);
// polynomial (FEM basis function) degree
int pp = 2;
const int nLocBas = pp + 1; // number of local basis, local means element
const int nFunc = pp * nElem + 1;
int IEN[nLocBas][nElem];
for (int ee = 0; ee < nElem; ee++)
{
for (int aa = 0; aa < nLocBas; aa++)
{
IEN[aa][ee] = ee * pp + aa;
}
}
const double hh = (omega_r - omega_l) / double(pp*nElem);
double x_coor[pp*nElem+1];
for (int ii = 0; ii < pp*nElem+1; ii++)
{
x_coor[ii] = omega_l + double(ii) * hh;
}
int ID[nFunc];
for (int ii = 0; ii < nFunc; ii++)
{
ID[ii] = ii;
}
ID[0] = -1; // assign the ID for the Dirichlet node to be -1
SparseMatrix <double,RowMajor> K(nFunc,nFunc);
VectorXd F(nFunc);
VectorXd u(nFunc);
for (int ee = 0; ee < nElem; ee++)
{
double k_ele[nLocBas][nLocBas];
double f_ele[nLocBas];
double x_ele[nLocBas];
for (int ii = 0; ii < nLocBas; ii++)
{
for (int jj = 0; jj < nLocBas; jj++)
{
k_ele[ii][jj] = 0.0;
}
f_ele[ii] = 0.0;
x_ele[ii] = 0.0;
}
for (int aa = 0; aa < nLocBas; aa++)
{
x_ele[aa] = x_coor[IEN[aa][ee]];
}
for (int qua = 0; qua < nqp; qua++)
{
double dx_dxi = 0.0;
double x_qua = 0.0;
for (int aa = 0; aa < nLocBas; aa++)
{
x_qua = x_qua + x_ele[aa] * PolyBasis(pp, aa+1, 0, qp[qua]);
dx_dxi = dx_dxi + x_ele[aa] * PolyBasis(pp, aa+1, 1, qp[qua]);
}
double dxi_dx = 1.0 / dx_dxi;
// element assembly for the K and F
for (int aa = 0; aa < nLocBas; aa++)
{
double Na = PolyBasis(pp, aa+1, 0, qp[qua]);
double Na_xi = PolyBasis(pp, aa+1, 1, qp[qua]);
f_ele[aa] = f_ele[aa] + wq[qua] * Func_source(x_qua) * Na * dx_dxi;
for (int bb = 0; bb < nLocBas; bb++)
{
double Nb_xi = PolyBasis(pp, bb+1, 1, qp[qua]);
k_ele[aa][bb] = k_ele[aa][bb] + wq[qua] * Na_xi * Nb_xi * dxi_dx;
}
}
}
for (int aa = 0; aa < nLocBas; aa++)
{
int AA = ID[IEN[aa][ee]];
if(AA >= 0)
{
F(AA) += f_ele[aa];
for (int bb = 0; bb < nLocBas; bb++)
{
int BB = IEN[bb][ee];
K.coeffRef(AA,BB) += k_ele[aa][bb];
}
}
else
{
K.coeffRef(IEN[aa][ee],IEN[aa][ee]) = 1.0;
F(IEN[aa][ee]) = g0;
}
}
}
F(nFunc-1) += gL;
SparseLU<SparseMatrix<double> > solver;
solver.compute(K);
u = solver.solve(F);
cout << "The number of elements = " << nElem << endl;
cout << "The degree of Legrange-shape function = " << pp << endl;
cout << "Displacement u = " << endl;
cout << u << endl;
}