# 引子

• 在不断细化的网格上的计算。数值计算通常要在不同的网格上进行，这样才能感受到精度。而且deal.II支持自适应网格，虽然这个例子中没有用到，但基础在这
• 读入非规则网格数据
• 计算优化
• debug模式，使用Assert宏
• 变系数Possion方程，使用预条件迭代求解器

a(x)u(x)u=1inΩ,=0onΩ

(aϕ,u)=(ϕ,1)ϕ

# 程序解析

 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  #include #include #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; 

step5类模板：

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19  template class Step5 { public: Step5 (); void run (); private: void setup_system (); void assemble_system (); void solve (); void output_results (const unsigned int cycle) 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  template class Coefficient : public Function { public: Coefficient () : Function() {} virtual double value (const Point &p, const unsigned int component = 0) const; virtual void value_list (const std::vector > &points, std::vector &values, const unsigned int component = 0) const; }; 

 1 2 3 4 5 6 7 8 9  template double Coefficient::value (const Point &p, const unsigned int / *component* /) const { if (p.square() < 0.5*0.5) return 20; else return 1; } 

 1 2 3 4 5 6 7  template void Coefficient::value_list (const std::vector > &points, std::vector &values, const unsigned int component) const { Assert (values.size() == points.size(), ExcDimensionMismatch (values.size(), points.size())); 

 1 2  Assert (component == 0, ExcIndexRange (component, 0, 1)); 

 1 2 3 4 5 6 7 8 9  const unsigned int n_points = points.size(); for (unsigned int i=0; i

 1 2 3 4 5  template Step5::Step5 () : fe (1), dof_handler (triangulation) {} 

 1 2 3 4 5 6 7 8 9 10 11 12 13 14  template void Step5::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 5 6 7 8 9 10 11 12  template void Step5::assemble_system () { QGauss quadrature_formula(2); FEValues fe_values (fe, quadrature_formula, update_values | update_gradients | update_quadrature_points | update_JxW_values); 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); 

 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  const Coefficient coefficient; std::vector coefficient_values (n_q_points); typename DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) { cell_matrix = 0; cell_rhs = 0; fe_values.reinit (cell); coefficient.value_list (fe_values.get_quadrature_points(), coefficient_values); for (unsigned int q_index=0; q_indexget_dof_indices (local_dof_indices); for (unsigned int i=0; i boundary_values; VectorTools::interpolate_boundary_values (dof_handler, 0, ZeroFunction(), boundary_values); MatrixTools::apply_boundary_values (boundary_values, system_matrix, solution, system_rhs); } 

 1 2 3 4 5 6 7 8 9 10 11 12 13  template void Step5::solve () { SolverControl solver_control (1000, 1e-12); SolverCG<> solver (solver_control); PreconditionSSOR<> preconditioner; preconditioner.initialize(system_matrix, 1.2); solver.solve (system_matrix, solution, system_rhs, preconditioner); 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 13 14 15 16 17 18 19  template void Step5::output_results (const unsigned int cycle) const { DataOut data_out; data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "solution"); data_out.build_patches (); DataOutBase::EpsFlags eps_flags; eps_flags.z_scaling = 4; eps_flags.azimut_angle = 40; eps_flags.turn_angle = 10; data_out.set_flags (eps_flags); std::ostringstream filename; filename << "solution-" << cycle << ".eps"; std::ofstream output (filename.str().c_str()); data_out.write_eps (output); } 

 1 2 3 4 5 6  template void Step5::run () { GridIn grid_in; grid_in.attach_triangulation (triangulation); std::ifstream input_file("circle-grid.inp"); 

run函数中直接读入网格文件，这里是个inp后缀的。因为该网格是二维的，所以Assert一下如果不是二维的，就抛出一个异常：

 1  Assert (dim==2, ExcInternalError()); 

 1  grid_in.read_ucd (input_file); 
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20  static const SphericalManifold boundary; triangulation.set_all_manifold_ids_on_boundary(0); triangulation.set_manifold (0, boundary); for (unsigned int cycle=0; cycle<6; ++cycle) { std::cout << "Cycle " << cycle << ':' << std::endl; if (cycle != 0) triangulation.refine_global (1); std::cout << " Number of active cells: " << triangulation.n_active_cells() << std::endl << " Total number of cells: " << triangulation.n_cells() << std::endl; setup_system (); assemble_system (); solve (); output_results (cycle); } } 

 1 2 3 4 5 6 7 8 9 10 11 12  int main () { Step5<2> laplace_problem_2d; laplace_problem_2d.run (); /* Coefficient<2> coefficient; std::vector > points (2); std::vector coefficient_values (1); coefficient.value_list (points, coefficient_values); */ return 0; } `

# 计算结果

• 本文已收录于以下专栏：

举报原因： 您举报文章：求解偏微分方程开源有限元软件deal.II学习--Step 5 色情 政治 抄袭 广告 招聘 骂人 其他 (最多只允许输入30个字)