📄 ex11.php
字号:
</pre></div><div class = "comment">The boundary values. <br><br>Set u = 1 on the top boundary, 0 everywhere else</div><div class ="fragment"><pre> const Real u_value = (yf > .99) ? 1. : 0.; </pre></div><div class = "comment">Set v = 0 everywhere</div><div class ="fragment"><pre> const Real v_value = 0.; </pre></div><div class = "comment">Find the node on the element matching this node onthe side. That defined where in the element matrixthe boundary condition will be applied.</div><div class ="fragment"><pre> for (unsigned int n=0; n<elem->n_nodes(); n++) if (elem->node(n) == side->node(ns)) {</pre></div><div class = "comment">Matrix contribution.</div><div class ="fragment"><pre> Kuu(n,n) += penalty; Kvv(n,n) += penalty; </pre></div><div class = "comment">Right-hand-side contribution.</div><div class ="fragment"><pre> Fu(n) += penalty*u_value; Fv(n) += penalty*v_value; } } // end face node loop } // end if (elem->neighbor(side) == NULL) } // end boundary condition section </pre></div><div class = "comment">The element matrix and right-hand-side are now builtfor this element. Add them to the global matrix andright-hand-side vector. The \p PetscMatrix::add_matrix()and \p PetscVector::add_vector() members do this for us.</div><div class ="fragment"><pre> system.matrix->add_matrix (Ke, dof_indices); system.rhs->add_vector (Fe, dof_indices); } // end of element loop </pre></div><div class = "comment">That's it.</div><div class ="fragment"><pre> return; }</pre></div><a name="nocomments"></a> <br><br><br> <h1> The program without comments: </h1> <pre> #include <iostream> #include <algorithm> #include <math.h> #include <B><FONT COLOR="#BC8F8F">"libmesh.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"mesh.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"mesh_generation.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"gmv_io.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"equation_systems.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"fe.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"quadrature_gauss.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"dof_map.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"sparse_matrix.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"numeric_vector.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"dense_matrix.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"dense_vector.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"linear_implicit_system.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"dense_submatrix.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"dense_subvector.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"elem.h"</FONT></B> <B><FONT COLOR="#228B22">void</FONT></B> assemble_stokes (EquationSystems& es, <B><FONT COLOR="#228B22">const</FONT></B> std::string& system_name); <B><FONT COLOR="#228B22">int</FONT></B> main (<B><FONT COLOR="#228B22">int</FONT></B> argc, <B><FONT COLOR="#228B22">char</FONT></B>** argv) { <B><FONT COLOR="#5F9EA0">libMesh</FONT></B>::init (argc, argv); { <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> dim = 2; Mesh mesh (dim); <B><FONT COLOR="#5F9EA0">MeshTools</FONT></B>::Generation::build_square (mesh, 15, 15, 0., 1., 0., 1., QUAD9); mesh.print_info(); EquationSystems equation_systems (mesh); { LinearImplicitSystem & system = equation_systems.add_system<LinearImplicitSystem> (<B><FONT COLOR="#BC8F8F">"Stokes"</FONT></B>); system.add_variable (<B><FONT COLOR="#BC8F8F">"u"</FONT></B>, SECOND); system.add_variable (<B><FONT COLOR="#BC8F8F">"v"</FONT></B>, SECOND); system.add_variable (<B><FONT COLOR="#BC8F8F">"p"</FONT></B>, FIRST); system.attach_assemble_function (assemble_stokes); equation_systems.init (); 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; equation_systems.print_info(); } equation_systems.get_system(<B><FONT COLOR="#BC8F8F">"Stokes"</FONT></B>).solve(); GMVIO(mesh).write_equation_systems (<B><FONT COLOR="#BC8F8F">"out.gmv"</FONT></B>, equation_systems); } <B><FONT COLOR="#A020F0">return</FONT></B> libMesh::close (); } <B><FONT COLOR="#228B22">void</FONT></B> assemble_stokes (EquationSystems& es, <B><FONT COLOR="#228B22">const</FONT></B> std::string& system_name) { assert (system_name == <B><FONT COLOR="#BC8F8F">"Stokes"</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(); LinearImplicitSystem & system = es.get_system<LinearImplicitSystem> (<B><FONT COLOR="#BC8F8F">"Stokes"</FONT></B>); <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> u_var = system.variable_number (<B><FONT COLOR="#BC8F8F">"u"</FONT></B>); <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> v_var = system.variable_number (<B><FONT COLOR="#BC8F8F">"v"</FONT></B>); <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> p_var = system.variable_number (<B><FONT COLOR="#BC8F8F">"p"</FONT></B>); FEType fe_vel_type = system.variable_type(u_var); FEType fe_pres_type = system.variable_type(p_var); AutoPtr<FEBase> fe_vel (FEBase::build(dim, fe_vel_type)); AutoPtr<FEBase> fe_pres (FEBase::build(dim, fe_pres_type)); QGauss qrule (dim, fe_vel_type.default_quadrature_order()); fe_vel->attach_quadrature_rule (&qrule); fe_pres->attach_quadrature_rule (&qrule); <B><FONT COLOR="#228B22">const</FONT></B> std::vector<Real>& JxW = fe_vel->get_JxW(); <B><FONT COLOR="#228B22">const</FONT></B> std::vector<std::vector<RealGradient> >& dphi = fe_vel->get_dphi(); <B><FONT COLOR="#228B22">const</FONT></B> std::vector<std::vector<Real> >& psi = fe_pres->get_phi(); <B><FONT COLOR="#228B22">const</FONT></B> DofMap & dof_map = system.get_dof_map(); DenseMatrix<Number> Ke; DenseVector<Number> Fe; DenseSubMatrix<Number> Kuu(Ke), Kuv(Ke), Kup(Ke), Kvu(Ke), Kvv(Ke), Kvp(Ke), Kpu(Ke), Kpv(Ke), Kpp(Ke); DenseSubVector<Number> Fu(Fe), Fv(Fe), Fp(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="#5F9EA0">std</FONT></B>::vector<<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B>> dof_indices_u; <B><FONT COLOR="#5F9EA0">std</FONT></B>::vector<<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B>> dof_indices_v; <B><FONT COLOR="#5F9EA0">std</FONT></B>::vector<<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B>> dof_indices_p; <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); dof_map.dof_indices (elem, dof_indices_u, u_var); dof_map.dof_indices (elem, dof_indices_v, v_var); dof_map.dof_indices (elem, dof_indices_p, p_var); <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> n_dofs = dof_indices.size(); <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> n_u_dofs = dof_indices_u.size(); <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> n_v_dofs = dof_indices_v.size(); <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> n_p_dofs = dof_indices_p.size(); fe_vel->reinit (elem); fe_pres->reinit (elem); Ke.resize (n_dofs, n_dofs); Fe.resize (n_dofs); Kuu.reposition (u_var*n_u_dofs, u_var*n_u_dofs, n_u_dofs, n_u_dofs); Kuv.reposition (u_var*n_u_dofs, v_var*n_u_dofs, n_u_dofs, n_v_dofs); Kup.reposition (u_var*n_u_dofs, p_var*n_u_dofs, n_u_dofs, n_p_dofs); Kvu.reposition (v_var*n_v_dofs, u_var*n_v_dofs, n_v_dofs, n_u_dofs); Kvv.reposition (v_var*n_v_dofs, v_var*n_v_dofs, n_v_dofs, n_v_dofs); Kvp.reposition (v_var*n_v_dofs, p_var*n_v_dofs, n_v_dofs, n_p_dofs); Kpu.reposition (p_var*n_u_dofs, u_var*n_u_dofs, n_p_dofs, n_u_dofs); Kpv.reposition (p_var*n_u_dofs, v_var*n_u_dofs, n_p_dofs, n_v_dofs); Kpp.reposition (p_var*n_u_dofs, p_var*n_u_dofs, n_p_dofs, n_p_dofs); Fu.reposition (u_var*n_u_dofs, n_u_dofs); Fv.reposition (v_var*n_u_dofs, n_v_dofs); Fp.reposition (p_var*n_u_dofs, n_p_dofs); <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++) { <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> i=0; i<n_u_dofs; 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<n_u_dofs; j++) Kuu(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][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<n_u_dofs; 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<n_p_dofs; j++) Kup(i,j) += -JxW[qp]*psi[j][qp]*dphi[i][qp](0); <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> i=0; i<n_v_dofs; 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<n_v_dofs; j++) Kvv(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][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<n_v_dofs; 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<n_p_dofs; j++) Kvp(i,j) += -JxW[qp]*psi[j][qp]*dphi[i][qp](1); <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> i=0; i<n_p_dofs; 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<n_u_dofs; j++) Kpu(i,j) += -JxW[qp]*psi[i][qp]*dphi[j][qp](0); <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> i=0; i<n_p_dofs; 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<n_v_dofs; j++) Kpv(i,j) += -JxW[qp]*psi[i][qp]*dphi[j][qp](1); } <I><FONT COLOR="#B22222">// end of the quadrature point qp-loop</FONT></I> { <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) { AutoPtr<Elem> side (elem->build_side(s)); <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> ns=0; ns<side->n_nodes(); ns++) { <B><FONT COLOR="#228B22">const</FONT></B> Real xf = side->point(ns)(0); <B><FONT COLOR="#228B22">const</FONT></B> Real yf = side->point(ns)(1); <B><FONT COLOR="#228B22">const</FONT></B> Real penalty = 1.e10; <B><FONT COLOR="#228B22">const</FONT></B> Real u_value = (yf > .99) ? 1. : 0.; <B><FONT COLOR="#228B22">const</FONT></B> Real v_value = 0.; <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> n=0; n<elem->n_nodes(); n++) <B><FONT COLOR="#A020F0">if</FONT></B> (elem->node(n) == side->node(ns)) { Kuu(n,n) += penalty; Kvv(n,n) += penalty; Fu(n) += penalty*u_value; Fv(n) += penalty*v_value; } } <I><FONT COLOR="#B22222">// end face node loop </FONT></I> } <I><FONT COLOR="#B22222">// end if (elem->neighbor(side) == NULL)</FONT></I> } <I><FONT COLOR="#B22222">// end boundary condition section </FONT></I> system.matrix->add_matrix (Ke, dof_indices); system.rhs->add_vector (Fe, dof_indices); } <I><FONT COLOR="#B22222">// end of element loop</FONT></I> <B><FONT COLOR="#A020F0">return</FONT></B>; }</pre> <a name="output"></a> <br><br><br> <h1> The console output of the program: </h1> <pre>**************************************************************** Running Example ./ex11-devel*************************************************************** Mesh Information: mesh_dimension()=2 spatial_dimension()=3 n_nodes()=961 n_elem()=225 n_local_elem()=225 n_active_elem()=225 n_subdomains()=1 n_processors()=1 processor_id()=0 EquationSystems n_systems()=1 System "Stokes" Type "LinearImplicit" Variables="u" "v" "p" Finite Element Types="LAGRANGE" "LAGRANGE" "LAGRANGE" Approximation Orders="SECOND" "SECOND" "FIRST" n_dofs()=2178 n_local_dofs()=2178 n_constrained_dofs()=0 n_vectors()=1 **************************************************************** Done Running Example ./ex11-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 + -