📄 ex3.php
字号:
points on the face.</div><div class ="fragment"><pre> const std::vector<Real>& JxW_face = fe_face->get_JxW(); </pre></div><div class = "comment">The XYZ locations (in physical space) of thequadrature points on the face. This is wherewe will interpolate the boundary value function.</div><div class ="fragment"><pre> const std::vector<Point >& qface_point = fe_face->get_xyz(); </pre></div><div class = "comment">Compute the shape function values on the elementface.</div><div class ="fragment"><pre> fe_face->reinit(elem, side); </pre></div><div class = "comment">Loop over the face quadrature points for integration.</div><div class ="fragment"><pre> for (unsigned int qp=0; qp<qface.n_points(); qp++) { </pre></div><div class = "comment">The location on the boundary of the currentface quadrature point.</div><div class ="fragment"><pre> const Real xf = qface_point[qp](0); const Real yf = qface_point[qp](1); </pre></div><div class = "comment">The penalty value. \frac{1}{\epsilon}in the discussion above.</div><div class ="fragment"><pre> const Real penalty = 1.e10; </pre></div><div class = "comment">The boundary value.</div><div class ="fragment"><pre> const Real value = exact_solution(xf, yf); </pre></div><div class = "comment">Matrix contribution of the L2 projection. </div><div class ="fragment"><pre> for (unsigned int i=0; i<phi_face.size(); i++) for (unsigned int j=0; j<phi_face.size(); j++) Ke(i,j) += JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp]; </pre></div><div class = "comment">Right-hand-side contribution of the L2projection.</div><div class ="fragment"><pre> for (unsigned int i=0; i<phi_face.size(); i++) Fe(i) += JxW_face[qp]*penalty*value*phi_face[i][qp]; } } } </pre></div><div class = "comment">We have now finished the quadrature point loop,and have therefore applied all the boundary conditions.<br><br>The element matrix and right-hand-side are now builtfor this element. Add them to the global matrix andright-hand-side vector. The SparseMatrix::add_matrix()and NumericVector::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); } </pre></div><div class = "comment">All done!</div><div class ="fragment"><pre> }</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">"linear_implicit_system.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">"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">"elem.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"dof_map.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"boundary_mesh.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"boundary_info.h"</FONT></B> <B><FONT COLOR="#228B22">void</FONT></B> assemble_poisson(EquationSystems& es, <B><FONT COLOR="#228B22">const</FONT></B> std::string& system_name); Real exact_solution (<B><FONT COLOR="#228B22">const</FONT></B> Real x, <B><FONT COLOR="#228B22">const</FONT></B> Real y, <B><FONT COLOR="#228B22">const</FONT></B> Real z = 0.); <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="#5F9EA0">std</FONT></B>::cout << <B><FONT COLOR="#BC8F8F">"Running "</FONT></B> << argv[0]; <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">int</FONT></B> i=1; i<argc; i++) <B><FONT COLOR="#5F9EA0">std</FONT></B>::cout << <B><FONT COLOR="#BC8F8F">" "</FONT></B> << argv[i]; <B><FONT COLOR="#5F9EA0">std</FONT></B>::cout << std::endl << std::endl; Mesh mesh (2); <B><FONT COLOR="#5F9EA0">MeshTools</FONT></B>::Generation::build_square (mesh, 15, 15, -1., 1., -1., 1., QUAD9); mesh.print_info(); EquationSystems equation_systems (mesh); equation_systems.add_system<LinearImplicitSystem> (<B><FONT COLOR="#BC8F8F">"Poisson"</FONT></B>); equation_systems.get_system(<B><FONT COLOR="#BC8F8F">"Poisson"</FONT></B>).add_variable(<B><FONT COLOR="#BC8F8F">"u"</FONT></B>, SECOND); equation_systems.get_system(<B><FONT COLOR="#BC8F8F">"Poisson"</FONT></B>).attach_assemble_function (assemble_poisson); equation_systems.init(); equation_systems.print_info(); equation_systems.get_system(<B><FONT COLOR="#BC8F8F">"Poisson"</FONT></B>).solve(); GMVIO (mesh).write_equation_systems (<B><FONT COLOR="#BC8F8F">"out.gmv"</FONT></B>, equation_systems); <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> i=0; i<3; i++) { here(); equation_systems.get_system(<B><FONT COLOR="#BC8F8F">"Poisson"</FONT></B>).solve(); BoundaryMesh boundary_mesh (mesh.mesh_dimension()-1); mesh.boundary_info->sync(boundary_mesh, false); } } <B><FONT COLOR="#A020F0">return</FONT></B> libMesh::close(); } <B><FONT COLOR="#228B22">void</FONT></B> assemble_poisson(EquationSystems& es, <B><FONT COLOR="#228B22">const</FONT></B> std::string& system_name) { assert (system_name == <B><FONT COLOR="#BC8F8F">"Poisson"</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">"Poisson"</FONT></B>); <B><FONT COLOR="#228B22">const</FONT></B> DofMap& dof_map = system.get_dof_map(); FEType fe_type = dof_map.variable_type(0); AutoPtr<FEBase> fe (FEBase::build(dim, fe_type)); QGauss qrule (dim, FIFTH); fe->attach_quadrature_rule (&qrule); AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type)); QGauss qface(dim-1, FIFTH); 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<Point>& q_point = fe->get_xyz(); <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<RealGradient> >& dphi = fe->get_dphi(); 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="#5F9EA0">MeshBase</FONT></B>::const_element_iterator el = mesh.elements_begin(); <B><FONT COLOR="#228B22">const</FONT></B> MeshBase::const_element_iterator end_el = mesh.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++) { <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++) <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]*(dphi[i][qp]*dphi[j][qp]); } { <B><FONT COLOR="#228B22">const</FONT></B> Real x = q_point[qp](0); <B><FONT COLOR="#228B22">const</FONT></B> Real y = q_point[qp](1); <B><FONT COLOR="#228B22">const</FONT></B> Real eps = 1.e-3; <B><FONT COLOR="#228B22">const</FONT></B> Real fxy = -(exact_solution(x,y-eps) + exact_solution(x,y+eps) + exact_solution(x-eps,y) + exact_solution(x+eps,y) - 4.*exact_solution(x,y))/eps/eps; <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]*fxy*phi[i][qp]; } } { <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> side=0; side<elem->n_sides(); side++) <B><FONT COLOR="#A020F0">if</FONT></B> (elem->neighbor(side) == NULL) { <B><FONT COLOR="#228B22">const</FONT></B> std::vector<std::vector<Real> >& phi_face = fe_face->get_phi(); <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<Point >& qface_point = fe_face->get_xyz(); fe_face->reinit(elem, side); <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> Real xf = qface_point[qp](0); <B><FONT COLOR="#228B22">const</FONT></B> Real yf = qface_point[qp](1); <B><FONT COLOR="#228B22">const</FONT></B> Real penalty = 1.e10; <B><FONT COLOR="#228B22">const</FONT></B> Real value = exact_solution(xf, yf); <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_face.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<phi_face.size(); j++) Ke(i,j) += JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[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<phi_face.size(); i++) Fe(i) += JxW_face[qp]*penalty*value*phi_face[i][qp]; } } } system.matrix->add_matrix (Ke, dof_indices); system.rhs->add_vector (Fe, dof_indices); } }</pre> <a name="output"></a> <br><br><br> <h1> The console output of the program: </h1> <pre>**************************************************************** Running Example ./ex3-devel*************************************************************** Running ./ex3-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 "Poisson" Type "LinearImplicit" Variables="u" Finite Element Types="LAGRANGE" Approximation Orders="SECOND" n_dofs()=961 n_local_dofs()=961 n_constrained_dofs()=0 n_vectors()=1[0] ex3.C, line 162, compiled Jun 6 2007 at 11:54:25[0] ex3.C, line 162, compiled Jun 6 2007 at 11:54:25[0] ex3.C, line 162, compiled Jun 6 2007 at 11:54:25 **************************************************************** Done Running Example ./ex3-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 + -