432人阅读 评论(0)

引子

deal.II有一个特性，叫作”维度无关的编程“。前面的例子中，很多类都有一个尖括号括起的数字的后缀。

程序解析

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24  #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include using namespace dealii; 

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20  template class Step4 { public: Step4 (); void run (); private: void make_grid (); void setup_system(); void assemble_system (); void solve (); void output_results () const; Triangulation triangulation; FE_Q fe; DoFHandler dof_handler; SparsityPattern sparsity_pattern; SparseMatrix system_matrix; Vector solution; Vector system_rhs; }; 

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16  template class RightHandSide : public Function { public: RightHandSide () : Function() {} virtual double value (const Point &p, const unsigned int component = 0) const; }; template class BoundaryValues : public Function { public: BoundaryValues () : Function() {} virtual double value (const Point &p, const unsigned int component = 0) const; }; 

 1 2 3 4 5 6 7 8 9  template double RightHandSide::value (const Point &p, const unsigned int / *component* /) const { double return_value = 0.0; for (unsigned int i=0; i

 1 2 3 4 5 6  template double BoundaryValues::value (const Point &p, const unsigned int / *component* /) const { return p.square(); } 

 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  template Step4::Step4 () : fe (1), dof_handler (triangulation) {} template void Step4::make_grid () { GridGenerator::hyper_cube (triangulation, -1, 1); triangulation.refine_global (4); std::cout << " Number of active cells: " << triangulation.n_active_cells() << std::endl << " Total number of cells: " << triangulation.n_cells() << std::endl; } template void Step4::setup_system () { dof_handler.distribute_dofs (fe); std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl; 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()); } 

 1 2 3 4  template void Step4::assemble_system () { QGauss quadrature_formula(2); 

 1  const RightHandSide right_hand_side; 

 1 2 3  FEValues fe_values (fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); 

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24  const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); Vector cell_rhs (dofs_per_cell); std::vector local_dof_indices (dofs_per_cell); typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) { fe_values.reinit (cell); cell_matrix = 0; cell_rhs = 0; for (unsigned int q_index=0; q_index

 1 2 3 4 5 6 7 8 9 10  cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; i

 1 2 3 4 5 6 7 8 9 10  std::map boundary_values; VectorTools::interpolate_boundary_values (dof_handler, 0, BoundaryValues(), boundary_values); MatrixTools::apply_boundary_values (boundary_values, system_matrix, solution, system_rhs); } 

 1 2 3 4 5 6 7 8 9 10 11  template void Step4::solve () { SolverControl solver_control (1000, 1e-12); SolverCG<> solver (solver_control); solver.solve (system_matrix, solution, system_rhs, PreconditionIdentity()); std::cout << " " << solver_control.last_step() << " CG iterations needed to obtain convergence." << std::endl; } 

 1 2 3 4 5 6 7 8 9 10 11 12  template void Step4::output_results () const { DataOut data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "solution"); data_out.build_patches (); std::ofstream output (dim == 2 ? "solution-2d.vtk" : "solution-3d.vtk"); data_out.write_vtk (output); } 

run函数无须多说：

 1 2 3 4 5 6 7 8 9 10  template void Step4::run () { std::cout << "Solving problem in " << dim << " space dimensions." << std::endl; make_grid(); setup_system (); assemble_system (); solve (); output_results (); } 

main函数如下：

 1 2 3 4 5 6 7 8 9 10 11 12 13  int main () { deallog.depth_console (0); { Step4<2> laplace_problem_2d; laplace_problem_2d.run (); } { Step4<3> laplace_problem_3d; laplace_problem_3d.run (); } return 0; } 

2d和3d的切换是如此从容，注意：这里用花括号将两者分开，是为了及时销毁变量和释放内存。

计算结果

 1 2 3 4 5 6 7 8 9 10  Solving problem in 2 space dimensions. Number of active cells: 256 Total number of cells: 341 Number of degrees of freedom: 289 26 CG iterations needed to obtain convergence. Solving problem in 3 space dimensions. Number of active cells: 4096 Total number of cells: 4681 Number of degrees of freedom: 4913 30 CG iterations needed to obtain convergence. `

2维计算结果为：

3维计算结果为：

0
0

* 以上用户言论只代表其个人观点，不代表CSDN网站的观点或立场
个人资料
• 访问：88742次
• 积分：1719
• 等级：
• 排名：千里之外
• 原创：34篇
• 转载：217篇
• 译文：19篇
• 评论：8条
GPU_CUDA异构并行计算
评论排行
阅读排行
最新评论
文史哲