算例3介绍了泊松方程在dealii中的求解
泊松方程为边界为0,右端项非零,可以转化为AU=F:
使用类Step3开始程序的运行:
class Step3
{
public:
Step3();
void run();
使用一些私有对象完成对应的功能:
private:
void make_grid();
void setup_system();
void assemble_system();
void solve();
void output_results() const;
跟step-1一样,定义网格节点变量:
Triangulation<2> triangulation;
FE_Q<2> fe;
DoFHandler<2> dof_handler;
由拉普拉斯方程离散化得到的系统矩阵的稀疏模式和值的变量:
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
方程右边的解和变量:
Vector<double> solution;
Vector<double> system_rhs;
传入类Step3中变量的值,默认构造函数会实现其他变量的值:
Step3::Step3()
: fe(1)
, dof_handler(triangulation)
{}
跟step-1一样,产生网格信息:
void Step3::make_grid()
{
GridGenerator::hyper_cube(triangulation, -1, 1);
triangulation.refine_global(5);
std::cout << "Number of active cells: " << triangulation.n_active_cells()
<< std::endl;
}
n_active_cells是指所有的网格,包括已加密的和原先的父单元。
然后将动态稀疏矩阵赋予sparsity_pattern中:
DynamicSparsityPattern dsp(dof_handler.n_dofs());
DoFTools::make_sparsity_pattern(dof_handler, dsp);
sparsity_pattern.copy_from(dsp);
为了区别对待稀疏矩阵与矩阵,转换稀疏矩阵为矩阵,这对于大规模计算非常重要:
system_matrix.reinit(sparsity_pattern);
最后是设置方程右手边的向量和解向量的大小:
solution.reinit(dof_handler.n_dofs());
system_rhs.reinit(dof_handler.n_dofs());
使用积分公式计算每个单元的积分:
QGauss<2> quadrature_formula(fe.degree + 1);
需要每个单元的不同的形函数、倒数信息、雅可比行列式的权重:
FEValues<2> fe_values(fe,
quadrature_formula,
update_values | update_gradients | update_JxW_values);
赋值一些常用的数值为变量,方便以后的重用以及更改,每个单元自由度个数以及求积分的点数目:
const unsigned int dofs_per_cell = fe.dofs_per_cell;
const unsigned int n_q_points = quadrature_formula.size();
使用矩阵计算每个单元的单纲,最后再放到总纲的稀疏矩阵中,以及右边的向量:
FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
Vector<double> cell_rhs(dofs_per_cell);
使用局部编号来计算自由度,当需要把局部单元自由度耦合到全局时,需要知道该单元的全局编号,使用types::global_dof_index来定义一个临时矩阵存储该信息:
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
使用循环来遍历全部单元,因为不修改cell,可以把变量cell设为const:
for (const auto &cell : dof_handler.active_cell_iterators())
{
因为每个单元的雅克比的导数是不一样的,所以在遍历时,需要重置变量cell,以及全局矩阵和右端项:
fe_values.reinit(cell);
cell_matrix = 0;
cell_rhs = 0;
开始遍历所有的积分点,由于每个积分点都与其他点有关,需要两层循环,使用i和j,使用shape_grid代表积分点处的形函数导数,JxW表示权重,对右端项做同样的操作:
for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
{
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
cell_matrix(i, j) +=
(fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
fe_values.JxW(q_index)); // dx
for (unsigned int i = 0; i < dofs_per_cell; ++i)
cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
1 * // f(x_q)
fe_values.JxW(q_index)); // dx
}
组总纲,需要单元自由度的总数:
cell->get_dof_indices(local_dof_indices);
遍历所有的单元,将单纲添加到总纲上:
for (unsigned int i = 0; i < dofs_per_cell; ++i)
for (unsigned int j = 0; j < dofs_per_cell; ++j)
system_matrix.add(local_dof_indices[i],
local_dof_indices[j],
cell_matrix(i, j));
右端项
for (unsigned int i = 0; i < dofs_per_cell; ++i)
system_rhs(local_dof_indices[i]) += cell_rhs(i);
}
需要把约束或者说边界条件加到方程中,否则方程有无数个解,首先获得边界的自由度和形函数,使用VectorTools::interpolate_boundary_values()插值边界函数,边界有很多种形式,不同的问题边界条件不一,需要单独添加,使用Functions::ZeroFunction(处处为0的方程)定义VectorTools::interpolate_boundary_values(),使用std::map绑定总纲和边界约束:
std::map<types::global_dof_index, double> boundary_values;
VectorTools::interpolate_boundary_values(dof_handler,
0,
Functions::ZeroFunction<2>(),
boundary_values);
MatrixTools::apply_boundary_values(boundary_values,
system_matrix,
solution,
system_rhs);
}
求解方程,收敛条件为运行1000次或者残差范数小于:10^12,使用SolverCG模板求解,参数使用默认,CG求解器有一个前置调节器作为它的第四个参数。我们还没有准备好深入研究这个问题,所以我们使用PreconditionIdentity作为预处理:
void Step3::solve()
{
SolverControl solver_control(1000, 1e-12);
SolverCG<Vector<double>> solver(solver_control);
solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
}
输出结果,把数据dof_handler和solution传递给data_out,把前端数据转化为中间数据格式(后端数据),:
void Step3::output_results() const
{
DataOut<2> data_out;
data_out.attach_dof_handler(dof_handler);
data_out.add_data_vector(solution, "solution");
data_out.build_patches();
#写入数据
std::ofstream output("solution.vtk");
data_out.write_vtk(output);
}
调用类Step3中的其他函数运行:
void Step3::run()
{
make_grid();
setup_system();
assemble_system();
solve();
output_results();
}
主函数:
int main()
{
deallog.depth_console(2);
Step3 laplace_problem;
laplace_problem.run();
return 0;
}