-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathboundary_value.cpp
More file actions
157 lines (146 loc) · 5.05 KB
/
Copy pathboundary_value.cpp
File metadata and controls
157 lines (146 loc) · 5.05 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
#include "ISOP2P1.h"
#include "preconditioner.h"
#include "functions.h"
#define DIM 2
void ISOP2P1::boundaryValueStokes(Vector<double> &x)
{
/// 各空间自由度.
unsigned int n_dof_v = fem_space_v.n_dof();
unsigned int n_dof_p = fem_space_p.n_dof();
unsigned int n_total_dof_v = 2 * n_dof_v;
const std::size_t * rowstart = sp_stokes.get_rowstart_indices();
const unsigned int * colnum = sp_stokes.get_column_numbers();
std::cout << "n_dof_v: " << n_dof_v << ", n_dof_p: " << n_dof_p << std::endl;
std::cout << "n_A: " << sp_stokes.n_rows() << ", m_A: " << sp_stokes.n_cols() << std::endl;
/// 遍历全部维度的速度节点.
for (unsigned int i = 0; i < n_total_dof_v; ++i)
{
/// 边界标志.
int bm = -1;
/// 判断一下是 x 方向还是 y 方向. 分别读取标志.
if (i < n_dof_v)
bm = fem_space_v.dofInfo(i).boundary_mark;
else
bm = fem_space_v.dofInfo(i - n_dof_v).boundary_mark;
if (bm == 0)
continue;
/// 对 Dirichelet 边界根据边界分别赋值. 注意同时还要区别 x 和
/// 方腔流边界条件.
if (bm == 1 || bm == 2 || bm == 4 || bm == 5)
x(i) = 0.0;
else if (bm == 3)
if (i < n_dof_v)
{
// /// 不包括顶端的两个端点,称为watertight cavity.
// x(i) = 1.0;
Regularized regular;
x(i) = regular.value(fem_space_v.dofInfo(i).interp_point);
}
else
{
x(i) = 0.0;
}
// /// poiseuille flow 边界条件设置
// if (bm == 1 || bm == 2 || bm == 4 || bm == 5)
// if (i < n_dof_v)
// {
// PoiseuilleVx poiseuille_vx(-1.0, 1.0);
// x(i) = poiseuille_vx.value(fem_space_v.dofInfo(i).interp_point);
// }
// else
// {
// PoiseuilleVy poiseuille_vy;
// x(i) = poiseuille_vy.value(fem_space_v.dofInfo(i - n_dof_v).interp_point);
// }
/// 右端项这样改, 如果该行和列其余元素均为零, 则在迭代中确
/// 保该数值解和边界一致.
// if (bm == 1 || bm == 2 || bm == 4 || bm == 5)
if (bm == 1 || bm == 2 || bm == 3 || bm == 4 || bm == 5)
{
rhs(i) = matrix.diag_element(i) * x(i);
/// 遍历 i 行.
for (unsigned int j = rowstart[i] + 1;
j < rowstart[i + 1]; ++j)
{
/// 第 j 个元素消成零(不是第 j 列!). 注意避开了对角元.
matrix.global_entry(j) -= matrix.global_entry(j);
/// 第 j 个元素是第 k 列.
unsigned int k = colnum[j];
/// 看看第 k 行的 i 列是否easymesh 为零元.
const unsigned int *p = std::find(&colnum[rowstart[k] + 1],
&colnum[rowstart[k + 1]],
i);
/// 如果是非零元. 则需要将这一项移动到右端项. 因为第 i 个未知量已知.
if (p != &colnum[rowstart[k + 1]])
{
/// 计算 k 行 i 列的存储位置.
unsigned int l = p - &colnum[rowstart[0]];
/// 移动到右端项. 等价于 r(k) = r(k) - x(i) * A(k, i).
rhs(k) -= matrix.global_entry(l)
* x(i);
/// 移完此项自然是零.
matrix.global_entry(l) -= matrix.global_entry(l);
}
}
}
}
std::cout << "boundary values for Stokes OK!" << std::endl;
};
void ISOP2P1::boundaryValueNS(Vector<double> &x)
{
/// 各空间自由度.
unsigned int n_dof_v = fem_space_v.n_dof();
unsigned int n_dof_p = fem_space_p.n_dof();
unsigned int n_total_dof_v = 2 * n_dof_v;
const std::size_t * rowstart = sp_stokes.get_rowstart_indices();
const unsigned int * colnum = sp_stokes.get_column_numbers();
for (unsigned int i = 0; i < n_total_dof_v; ++i)
{
/// 边界标志.
int bm = -1;
/// 判断一下是 x 方向还是 y 方向. 分别读取标志.
if (i < n_dof_v)
bm = fem_space_v.dofInfo(i).boundary_mark;
else
bm = fem_space_v.dofInfo(i - n_dof_v).boundary_mark;
if (bm == 0)
continue;
// /// 对全部 Dirichlet 边界.
if (bm == 2 || bm == 3 || bm == 5 || bm == 1 || bm == 4 || bm == 11)
// if (bm < 10 && bm > 0 && bm != 6)
// if (bm == 1 || bm == 2 || bm == 3 || bm == 4)
{
/// 数值解对应点按成驱动速度.
x(i) = 0.0;
/// 右端项这样改, 如果该行和列其余元素均为零, 则在迭代中确
/// 保该数值解和边界一致.
rhs(i) = matrix.diag_element(i) * x(i);
/// 遍历 i 行.
for (unsigned int j = rowstart[i] + 1;
j < rowstart[i + 1]; ++j)
{
/// 第 j 个元素消成零(不是第 j 列!). 注意避开了对角元.
matrix.global_entry(j) -= matrix.global_entry(j);
/// 第 j 个元素是第 k 列.
unsigned int k = colnum[j];
/// 看看第 k 行的 i 列是否为零元.
const unsigned int *p = std::find(&colnum[rowstart[k] + 1],
&colnum[rowstart[k + 1]],
i);
/// 如果是非零元. 则需要将这一项移动到右端项. 因为第 i 个未知量已知.
if (p != &colnum[rowstart[k + 1]])
{
/// 计算 k 行 i 列的存储位置.
unsigned int l = p - &colnum[rowstart[0]];
/// 移动到右端项. 等价于 r(k) = r(k) - x(i) * A(k, i).
rhs(k) -= matrix.global_entry(l)
* x(i);
/// 移完此项自然是零.
matrix.global_entry(l) -= matrix.global_entry(l);
}
}
}
}
std::cout << "boundary apply to NS OK!" << std::endl;
};
#undef DIM