📄 ex10.php
字号:
<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> init_timestep = 0; <B><FONT COLOR="#A020F0">if</FONT></B>(command_line.search(<B><FONT COLOR="#BC8F8F">"-init_timestep"</FONT></B>)) init_timestep = command_line.next(0); <B><FONT COLOR="#A020F0">else</FONT></B> { <B><FONT COLOR="#5F9EA0">std</FONT></B>::cerr << <B><FONT COLOR="#BC8F8F">"ERROR: Initial timestep not specified\n"</FONT></B> << std::endl; error(); } <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> n_timesteps = 0; <B><FONT COLOR="#A020F0">if</FONT></B>(command_line.search(<B><FONT COLOR="#BC8F8F">"-n_timesteps"</FONT></B>)) n_timesteps = command_line.next(0); <B><FONT COLOR="#A020F0">else</FONT></B> { <B><FONT COLOR="#5F9EA0">std</FONT></B>::cerr << <B><FONT COLOR="#BC8F8F">"ERROR: Number of timesteps not specified\n"</FONT></B> << std::endl; error(); } Mesh mesh (2); EquationSystems equation_systems (mesh); MeshRefinement mesh_refinement (mesh); <B><FONT COLOR="#A020F0">if</FONT></B>(!read_solution) { mesh.read (<B><FONT COLOR="#BC8F8F">"mesh.xda"</FONT></B>); <B><FONT COLOR="#A020F0">if</FONT></B>(!read_solution) mesh_refinement.uniformly_refine (5); mesh.print_info(); TransientLinearImplicitSystem & system = equation_systems.add_system<TransientLinearImplicitSystem> (<B><FONT COLOR="#BC8F8F">"Convection-Diffusion"</FONT></B>); system.add_variable (<B><FONT COLOR="#BC8F8F">"u"</FONT></B>, FIRST); system.attach_assemble_function (assemble_cd); system.attach_init_function (init_cd); equation_systems.init (); } <B><FONT COLOR="#A020F0">else</FONT></B> { mesh.read(<B><FONT COLOR="#BC8F8F">"saved_mesh.xda"</FONT></B>); mesh.print_info(); equation_systems.read(<B><FONT COLOR="#BC8F8F">"saved_solution.xda"</FONT></B>, libMeshEnums::READ); TransientLinearImplicitSystem & system = equation_systems.get_system<TransientLinearImplicitSystem> (<B><FONT COLOR="#BC8F8F">"Convection-Diffusion"</FONT></B>); system.update(); system.attach_assemble_function (assemble_cd); } equation_systems.print_info(); equation_systems.parameters.set<<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B>> (<B><FONT COLOR="#BC8F8F">"linear solver maximum iterations"</FONT></B>) = 250; equation_systems.parameters.set<Real> (<B><FONT COLOR="#BC8F8F">"linear solver tolerance"</FONT></B>) = TOLERANCE; <B><FONT COLOR="#A020F0">if</FONT></B>(!read_solution) GMVIO(mesh).write_equation_systems (<B><FONT COLOR="#BC8F8F">"out.gmv.000"</FONT></B>, equation_systems); <B><FONT COLOR="#A020F0">else</FONT></B> GMVIO(mesh).write_equation_systems (<B><FONT COLOR="#BC8F8F">"solution_read_in.gmv"</FONT></B>, equation_systems); equation_systems.parameters.set<RealVectorValue>(<B><FONT COLOR="#BC8F8F">"velocity"</FONT></B>) = RealVectorValue (0.8, 0.8); <B><FONT COLOR="#228B22">const</FONT></B> Real dt = 0.025; Real time = init_timestep*dt; <B><FONT COLOR="#A020F0">for</FONT></B>(<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> t_step=init_timestep; t_step<(init_timestep+n_timesteps); t_step++) { time += dt; equation_systems.parameters.set<Real> (<B><FONT COLOR="#BC8F8F">"time"</FONT></B>) = time; equation_systems.parameters.set<Real> (<B><FONT COLOR="#BC8F8F">"dt"</FONT></B>) = dt; <B><FONT COLOR="#5F9EA0">std</FONT></B>::cout << <B><FONT COLOR="#BC8F8F">" Solving time step "</FONT></B>; { OStringStream out; OSSInt(out,2,t_step); out << <B><FONT COLOR="#BC8F8F">", time="</FONT></B>; OSSRealzeroleft(out,6,3,time); out << <B><FONT COLOR="#BC8F8F">"..."</FONT></B>; <B><FONT COLOR="#5F9EA0">std</FONT></B>::cout << out.str() << std::endl; } TransientLinearImplicitSystem & system = equation_systems.get_system<TransientLinearImplicitSystem>(<B><FONT COLOR="#BC8F8F">"Convection-Diffusion"</FONT></B>); *system.old_local_solution = *system.current_local_solution; <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> max_r_steps = 2; <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> r_step=0; r_step<max_r_steps; r_step++) { system.solve(); <B><FONT COLOR="#A020F0">if</FONT></B> (r_step+1 != max_r_steps) { <B><FONT COLOR="#5F9EA0">std</FONT></B>::cout << <B><FONT COLOR="#BC8F8F">" Refining the mesh..."</FONT></B> << std::endl; ErrorVector error; KellyErrorEstimator error_estimator; error_estimator.estimate_error (system, error); mesh_refinement.refine_fraction() = 0.80; mesh_refinement.coarsen_fraction() = 0.07; mesh_refinement.max_h_level() = 5; mesh_refinement.flag_elements_by_error_fraction (error); mesh_refinement.refine_and_coarsen_elements(); equation_systems.reinit (); } } <B><FONT COLOR="#A020F0">if</FONT></B> ( (t_step+1)%10 == 0) { OStringStream file_name; file_name << <B><FONT COLOR="#BC8F8F">"out.gmv."</FONT></B>; OSSRealzeroright(file_name,3,0,t_step+1); GMVIO(mesh).write_equation_systems (file_name.str(), equation_systems); } } <B><FONT COLOR="#A020F0">if</FONT></B>(!read_solution) { mesh.write(<B><FONT COLOR="#BC8F8F">"saved_mesh.xda"</FONT></B>); equation_systems.write(<B><FONT COLOR="#BC8F8F">"saved_solution.xda"</FONT></B>, libMeshEnums::WRITE); } } #endif <I><FONT COLOR="#B22222">// #ifndef ENABLE_AMR</FONT></I> <B><FONT COLOR="#A020F0">return</FONT></B> libMesh::close (); } <B><FONT COLOR="#228B22">void</FONT></B> init_cd (EquationSystems& es, <B><FONT COLOR="#228B22">const</FONT></B> std::string& system_name) { assert (system_name == <B><FONT COLOR="#BC8F8F">"Convection-Diffusion"</FONT></B>); TransientLinearImplicitSystem & system = es.get_system<TransientLinearImplicitSystem>(<B><FONT COLOR="#BC8F8F">"Convection-Diffusion"</FONT></B>); es.parameters.set<Real> (<B><FONT COLOR="#BC8F8F">"time"</FONT></B>) = 0; system.project_solution(exact_value, NULL, es.parameters); } <B><FONT COLOR="#228B22">void</FONT></B> assemble_cd (EquationSystems& es, <B><FONT COLOR="#228B22">const</FONT></B> std::string& system_name) { #ifdef ENABLE_AMR assert (system_name == <B><FONT COLOR="#BC8F8F">"Convection-Diffusion"</FONT></B>); <B><FONT COLOR="#228B22">const</FONT></B> Mesh& mesh = es.get_mesh(); <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> dim = mesh.mesh_dimension(); TransientLinearImplicitSystem & system = es.get_system<TransientLinearImplicitSystem> (<B><FONT COLOR="#BC8F8F">"Convection-Diffusion"</FONT></B>); FEType fe_type = system.variable_type(0); AutoPtr<FEBase> fe (FEBase::build(dim, fe_type)); AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type)); QGauss qrule (dim, fe_type.default_quadrature_order()); QGauss qface (dim-1, fe_type.default_quadrature_order()); fe->attach_quadrature_rule (&qrule); fe_face->attach_quadrature_rule (&qface); <B><FONT COLOR="#228B22">const</FONT></B> std::vector<Real>& JxW = fe->get_JxW(); <B><FONT COLOR="#228B22">const</FONT></B> std::vector<Real>& JxW_face = fe_face->get_JxW(); <B><FONT COLOR="#228B22">const</FONT></B> std::vector<std::vector<Real> >& phi = fe->get_phi(); <B><FONT COLOR="#228B22">const</FONT></B> std::vector<std::vector<Real> >& psi = fe_face->get_phi(); <B><FONT COLOR="#228B22">const</FONT></B> std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi(); <B><FONT COLOR="#228B22">const</FONT></B> std::vector<Point>& qface_points = fe_face->get_xyz(); <B><FONT COLOR="#228B22">const</FONT></B> DofMap& dof_map = system.get_dof_map(); DenseMatrix<Number> Ke; DenseVector<Number> Fe; <B><FONT COLOR="#5F9EA0">std</FONT></B>::vector<<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B>> dof_indices; <B><FONT COLOR="#228B22">const</FONT></B> RealVectorValue velocity = es.parameters.get<RealVectorValue> (<B><FONT COLOR="#BC8F8F">"velocity"</FONT></B>); <B><FONT COLOR="#228B22">const</FONT></B> Real dt = es.parameters.get<Real> (<B><FONT COLOR="#BC8F8F">"dt"</FONT></B>); <B><FONT COLOR="#228B22">const</FONT></B> Real time = es.parameters.get<Real> (<B><FONT COLOR="#BC8F8F">"time"</FONT></B>); <B><FONT COLOR="#5F9EA0">MeshBase</FONT></B>::const_element_iterator el = mesh.active_local_elements_begin(); <B><FONT COLOR="#228B22">const</FONT></B> MeshBase::const_element_iterator end_el = mesh.active_local_elements_end(); <B><FONT COLOR="#A020F0">for</FONT></B> ( ; el != end_el; ++el) { <B><FONT COLOR="#228B22">const</FONT></B> Elem* elem = *el; dof_map.dof_indices (elem, dof_indices); fe->reinit (elem); Ke.resize (dof_indices.size(), dof_indices.size()); Fe.resize (dof_indices.size()); <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> qp=0; qp<qrule.n_points(); qp++) { Number u_old = 0.; Gradient grad_u_old; <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> l=0; l<phi.size(); l++) { u_old += phi[l][qp]*system.old_solution (dof_indices[l]); grad_u_old.add_scaled (dphi[l][qp],system.old_solution (dof_indices[l])); } <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> i=0; i<phi.size(); i++) { Fe(i) += JxW[qp]*( u_old*phi[i][qp] + -.5*dt*( (grad_u_old*velocity)*phi[i][qp] + 0.01*(grad_u_old*dphi[i][qp])) ); <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> j=0; j<phi.size(); j++) { Ke(i,j) += JxW[qp]*( phi[i][qp]*phi[j][qp] + .5*dt*( (velocity*dphi[j][qp])*phi[i][qp] + 0.01*(dphi[i][qp]*dphi[j][qp])) ); } } } { <B><FONT COLOR="#228B22">const</FONT></B> Real penalty = 1.e10; <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> s=0; s<elem->n_sides(); s++) <B><FONT COLOR="#A020F0">if</FONT></B> (elem->neighbor(s) == NULL) { fe_face->reinit(elem,s); <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> qp=0; qp<qface.n_points(); qp++) { <B><FONT COLOR="#228B22">const</FONT></B> Number value = exact_solution (qface_points[qp](0), qface_points[qp](1), time); <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> i=0; i<psi.size(); i++) Fe(i) += penalty*JxW_face[qp]*value*psi[i][qp]; <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> i=0; i<psi.size(); i++) <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> j=0; j<psi.size(); j++) Ke(i,j) += penalty*JxW_face[qp]*psi[i][qp]*psi[j][qp]; } } } dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices); system.matrix->add_matrix (Ke, dof_indices); system.rhs->add_vector (Fe, dof_indices); } #endif <I><FONT COLOR="#B22222">// #ifdef ENABLE_AMR</FONT></I> }</pre> <a name="output"></a> <br><br><br> <h1> The console output of the program: </h1> <pre>**************************************************************** Running Example ./ex10-devel*************************************************************** Usage: ./ex10-devel -init_timestep 0OR ./ex10-devel -read_solution -init_timestep 26Running: ./ex10-devel -n_timesteps 25 -init_timestep 0 Mesh Information: mesh_dimension()=2 spatial_dimension()=3 n_nodes()=6273 n_elem()=13650 n_local_elem()=13650 n_active_elem()=10240 n_subdomains()=1 n_processors()=1 processor_id()=0 EquationSystems n_systems()=1 System "Convection-Diffusion" Type "TransientLinearImplicit" Variables="u" Finite Element Types="LAGRANGE" Approximation Orders="FIRST" n_dofs()=6273 n_local_dofs()=6273 n_constrained_dofs()=0 n_vectors()=3 Solving time step 0, time=0.0250... Refining the mesh... Solving time step 1, time=0.0500... Refining the mesh... Solving time step 2, time=0.0750... Refining the mesh... Solving time step 3, time=0.1000... Refining the mesh... Solving time step 4, time=0.1250... Refining the mesh... Solving time step 5, time=0.1500... Refining the mesh... Solving time step 6, time=0.1750... Refining the mesh... Solving time step 7, time=0.2000... Refining the mesh... Solving time step 8, time=0.2250... Refining the mesh... Solving time step 9, time=0.2500... Refining the mesh... Solving time step 10, time=0.2750... Refining the mesh... Solving time step 11, time=0.3000... Refining the mesh... Solving time step 12, time=0.3250... Refining the mesh... Solving time step 13, time=0.3500... Refining the mesh... Solving time step 14, time=0.3750... Refining the mesh... Solving time step 15, time=0.4000... Refining the mesh... Solving time step 16, time=0.4250... Refining the mesh... Solving time step 17, time=0.4500... Refining the mesh... Solving time step 18, time=0.4750... Refining the mesh... Solving time step 19, time=0.5000... Refining the mesh... Solving time step 20, time=0.5250... Refining the mesh... Solving time step 21, time=0.5500... Refining the mesh... Solving time step 22, time=0.5750... Refining the mesh... Solving time step 23, time=0.6000... Refining the mesh... Solving time step 24, time=0.6250... Refining the mesh... ***** Finished first 25 steps, now read in saved solution and continue ***** Usage: ./ex10-devel -init_timestep 0OR ./ex10-devel -read_solution -init_timestep 26Running: ./ex10-devel -read_solution -n_timesteps 25 -init_timestep 25 Mesh Information: mesh_dimension()=2 spatial_dimension()=3 n_nodes()=713 n_elem()=1018 n_local_elem()=1018 n_active_elem()=766 n_subdomains()=1 n_processors()=1 processor_id()=0 EquationSystems n_systems()=1 System "Convection-Diffusion" Type "TransientLinearImplicit" Variables="u" Finite Element Types="LAGRANGE" Approximation Orders="FIRST" n_dofs()=713 n_local_dofs()=713 n_constrained_dofs()=122 n_vectors()=3 Solving time step 25, time=0.6500... Refining the mesh... Solving time step 26, time=0.6750... Refining the mesh... Solving time step 27, time=0.7000... Refining the mesh... Solving time step 28, time=0.7250... Refining the mesh... Solving time step 29, time=0.7500... Refining the mesh... Solving time step 30, time=0.7750... Refining the mesh... Solving time step 31, time=0.8000... Refining the mesh... Solving time step 32, time=0.8250... Refining the mesh... Solving time step 33, time=0.8500... Refining the mesh... Solving time step 34, time=0.8750... Refining the mesh... Solving time step 35, time=0.9000... Refining the mesh... Solving time step 36, time=0.9250... Refining the mesh... Solving time step 37, time=0.9500... Refining the mesh... Solving time step 38, time=0.9750... Refining the mesh... Solving time step 39, time=1.0000... Refining the mesh... Solving time step 40, time=1.0300... Refining the mesh... Solving time step 41, time=1.0500... Refining the mesh... Solving time step 42, time=1.0700... Refining the mesh... Solving time step 43, time=1.1000... Refining the mesh... Solving time step 44, time=1.1200... Refining the mesh... Solving time step 45, time=1.1500... Refining the mesh... Solving time step 46, time=1.1700... Refining the mesh... Solving time step 47, time=1.2000... Refining the mesh... Solving time step 48, time=1.2200... Refining the mesh... Solving time step 49, time=1.2500... Refining the mesh... **************************************************************** Done Running Example ./ex10-devel***************************************************************</pre></div><?php make_footer() ?></body></html><?php if (0) { ?>\#Local Variables:\#mode: html\#End:<?php } ?>⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -