📄 ex6.php
字号:
<div class ="fragment"><pre> for (unsigned int i=0; i<n_sf; i++) for (unsigned int j=0; j<n_sf; j++) {</pre></div><div class = "comment">(ndt*Ht + nHt*d) * nH </div><div class ="fragment"><pre> Ke(i,j) += ( // ( ( // ( dweight[qp] * phi[i][qp] // Point * Real = Point + // + dphi[i][qp] * weight[qp] // Point * Real = Point ) * dphi[j][qp] // ) * Point = Real ) * JxW[qp]; // ) * Real = Real </pre></div><div class = "comment">(d*Ht*nmut*nH - ndt*nmu*Ht*H - d*nHt*nmu*H)</div><div class ="fragment"><pre> Ce(i,j) += ( // ( (dphase[qp] * dphi[j][qp]) // (Point * Point) = Real * weight[qp] * phi[i][qp] // * Real * Real = Real - // - (dweight[qp] * dphase[qp]) // (Point * Point) = Real * phi[i][qp] * phi[j][qp] // * Real * Real = Real - // - (dphi[i][qp] * dphase[qp]) // (Point * Point) = Real * weight[qp] * phi[j][qp] // * Real * Real = Real ) * JxW[qp]; // ) * Real = Real </pre></div><div class = "comment">(d*Ht*H * (1 - nmut*nmu))</div><div class ="fragment"><pre> Me(i,j) += ( // ( (1. - (dphase[qp] * dphase[qp])) // (Real - (Point * Point)) = Real * phi[i][qp] * phi[j][qp] * weight[qp] // * Real * Real * Real = Real ) * JxW[qp]; // ) * Real = Real } // end of the matrix summation loop } // end of quadrature point loop </pre></div><div class = "comment">The element matrices are now built for this element. Collect them in Ke, and then add them to the global matrix. The \p SparseMatrix::add_matrix() member does this for us.</div><div class ="fragment"><pre> Ke.add(1./speed , Ce); Ke.add(1./(speed*speed), Me); system.matrix->add_matrix (Ke, dof_indices); } // end of element loop </pre></div><div class = "comment">Note that we have not applied any boundary conditions so far.Here we apply a unit load at the node located at (0,0,0).</div><div class ="fragment"><pre> {</pre></div><div class = "comment">Number of nodes in the mesh. </div><div class ="fragment"><pre> const unsigned int n_nodes = mesh.n_nodes(); for (unsigned int n_cnt=0; n_cnt<n_nodes; n_cnt++) { </pre></div><div class = "comment">Get a reference to the current node.</div><div class ="fragment"><pre> const Node& curr_node = mesh.node(n_cnt); </pre></div><div class = "comment">Check the location of the current node.</div><div class ="fragment"><pre> if (fabs(curr_node(0)) < TOLERANCE && fabs(curr_node(1)) < TOLERANCE && fabs(curr_node(2)) < TOLERANCE) {</pre></div><div class = "comment">The global number of the respective degree of freedom.</div><div class ="fragment"><pre> unsigned int dn = curr_node.dof_number(0,0,0); system.rhs->add (dn, 1.); } } } #else </pre></div><div class = "comment">dummy assert </div><div class ="fragment"><pre> assert(es.get_mesh().mesh_dimension() != 1); #endif //ifdef ENABLE_INFINITE_ELEMENTS </pre></div><div class = "comment">All done! </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">"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">"inf_fe.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"inf_elem_builder.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">"dof_map.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"node.h"</FONT></B> #include <B><FONT COLOR="#BC8F8F">"elem.h"</FONT></B> <B><FONT COLOR="#228B22">void</FONT></B> assemble_wave (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); #ifndef ENABLE_INFINITE_ELEMENTS <B><FONT COLOR="#5F9EA0">std</FONT></B>::cerr << <B><FONT COLOR="#BC8F8F">"ERROR: This example requires the library to be compiled with Infinite Element support!"</FONT></B> << std::endl; here(); <B><FONT COLOR="#A020F0">return</FONT></B> 0; #<B><FONT COLOR="#A020F0">else</FONT></B> { <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> dim = 3; <B><FONT COLOR="#5F9EA0">std</FONT></B>::cout << <B><FONT COLOR="#BC8F8F">"Running ex6 with dim = "</FONT></B> << dim << std::endl << std::endl; Mesh mesh (dim); <B><FONT COLOR="#5F9EA0">MeshTools</FONT></B>::Generation::build_cube (mesh, 4, 4, 4, -1., 1., -1., 1., -1., 1., HEX8); mesh.print_info(); GMVIO(mesh).write (<B><FONT COLOR="#BC8F8F">"orig_mesh.gmv"</FONT></B>); InfElemBuilder builder(mesh); builder.build_inf_elem(true); mesh.print_info(); GMVIO(mesh).write (<B><FONT COLOR="#BC8F8F">"ifems_added.gmv"</FONT></B>); mesh.find_neighbors(); EquationSystems equation_systems (mesh); { equation_systems.add_system<LinearImplicitSystem> (<B><FONT COLOR="#BC8F8F">"Wave"</FONT></B>); FEType fe_type(FIRST); equation_systems.get_system(<B><FONT COLOR="#BC8F8F">"Wave"</FONT></B>).add_variable(<B><FONT COLOR="#BC8F8F">"p"</FONT></B>, fe_type); equation_systems.get_system(<B><FONT COLOR="#BC8F8F">"Wave"</FONT></B>).attach_assemble_function (assemble_wave); equation_systems.parameters.set<Real>(<B><FONT COLOR="#BC8F8F">"speed"</FONT></B>) = 1.; equation_systems.parameters.set<Real>(<B><FONT COLOR="#BC8F8F">"fluid density"</FONT></B>) = 1.; equation_systems.init(); equation_systems.print_info(); } equation_systems.get_system(<B><FONT COLOR="#BC8F8F">"Wave"</FONT></B>).solve(); equation_systems.write (<B><FONT COLOR="#BC8F8F">"eqn_sys.dat"</FONT></B>, libMeshEnums::WRITE); } <B><FONT COLOR="#A020F0">return</FONT></B> libMesh::close (); #endif <I><FONT COLOR="#B22222">// else part of ifndef ENABLE_INFINITE_ELEMENTS</FONT></I> } <B><FONT COLOR="#228B22">void</FONT></B> assemble_wave(EquationSystems& es, <B><FONT COLOR="#228B22">const</FONT></B> std::string& system_name) { assert (system_name == <B><FONT COLOR="#BC8F8F">"Wave"</FONT></B>); #ifdef ENABLE_INFINITE_ELEMENTS <B><FONT COLOR="#228B22">const</FONT></B> Mesh& mesh = es.get_mesh(); LinearImplicitSystem & system = es.get_system<LinearImplicitSystem>(<B><FONT COLOR="#BC8F8F">"Wave"</FONT></B>); <B><FONT COLOR="#228B22">const</FONT></B> DofMap& dof_map = system.get_dof_map(); <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(); <B><FONT COLOR="#228B22">const</FONT></B> Real speed = es.parameters.get<Real>(<B><FONT COLOR="#BC8F8F">"speed"</FONT></B>); <B><FONT COLOR="#228B22">const</FONT></B> FEType& fe_type = dof_map.variable_type(0); AutoPtr<FEBase> fe (FEBase::build(dim, fe_type)); AutoPtr<FEBase> inf_fe (FEBase::build_InfFE(dim, fe_type)); QGauss qrule (dim, SECOND); fe->attach_quadrature_rule (&qrule); inf_fe->attach_quadrature_rule (&qrule); DenseMatrix<Number> Ke; DenseMatrix<Number> Ce; DenseMatrix<Number> Me; 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.local_elements_begin(); <B><FONT COLOR="#228B22">const</FONT></B> MeshBase::const_element_iterator end_el = mesh.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); FEBase* cfe=NULL; <B><FONT COLOR="#A020F0">if</FONT></B> (elem->infinite()) { cfe = inf_fe.get(); } <B><FONT COLOR="#A020F0">else</FONT></B> { cfe = fe.get(); { Fe.resize (dof_indices.size()); system.rhs->add_vector (Fe, dof_indices); } <I><FONT COLOR="#B22222">// end boundary condition section </FONT></I> } <I><FONT COLOR="#B22222">// else ( if (elem->infinite())) )</FONT></I> <B><FONT COLOR="#228B22">const</FONT></B> std::vector<Real>& JxW = cfe->get_JxW(); <B><FONT COLOR="#228B22">const</FONT></B> std::vector<std::vector<Real> >& phi = cfe->get_phi(); <B><FONT COLOR="#228B22">const</FONT></B> std::vector<std::vector<RealGradient> >& dphi = cfe->get_dphi(); <B><FONT COLOR="#228B22">const</FONT></B> std::vector<RealGradient>& dphase = cfe->get_dphase(); <B><FONT COLOR="#228B22">const</FONT></B> std::vector<Real>& weight = cfe->get_Sobolev_weight(); <B><FONT COLOR="#228B22">const</FONT></B> std::vector<RealGradient>& dweight = cfe->get_Sobolev_dweight(); cfe->reinit (elem); Ke.resize (dof_indices.size(), dof_indices.size()); Ce.resize (dof_indices.size(), dof_indices.size()); Me.resize (dof_indices.size(), dof_indices.size()); <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> max_qp = cfe->n_quadrature_points(); <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> qp=0; qp<max_qp; qp++) { <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> n_sf = cfe->n_shape_functions(); <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_sf; 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_sf; j++) { Ke(i,j) += ( <I><FONT COLOR="#B22222">// ( </FONT></I> ( <I><FONT COLOR="#B22222">// ( </FONT></I> dweight[qp] * phi[i][qp] <I><FONT COLOR="#B22222">// Point * Real = Point </FONT></I> + <I><FONT COLOR="#B22222">// + </FONT></I> dphi[i][qp] * weight[qp] <I><FONT COLOR="#B22222">// Point * Real = Point </FONT></I> ) * dphi[j][qp] <I><FONT COLOR="#B22222">// ) * Point = Real </FONT></I> ) * JxW[qp]; <I><FONT COLOR="#B22222">// ) * Real = Real </FONT></I> Ce(i,j) += ( <I><FONT COLOR="#B22222">// ( </FONT></I> (dphase[qp] * dphi[j][qp]) <I><FONT COLOR="#B22222">// (Point * Point) = Real </FONT></I> * weight[qp] * phi[i][qp] <I><FONT COLOR="#B22222">// * Real * Real = Real </FONT></I> - <I><FONT COLOR="#B22222">// - </FONT></I> (dweight[qp] * dphase[qp]) <I><FONT COLOR="#B22222">// (Point * Point) = Real </FONT></I> * phi[i][qp] * phi[j][qp] <I><FONT COLOR="#B22222">// * Real * Real = Real </FONT></I> - <I><FONT COLOR="#B22222">// - </FONT></I> (dphi[i][qp] * dphase[qp]) <I><FONT COLOR="#B22222">// (Point * Point) = Real </FONT></I> * weight[qp] * phi[j][qp] <I><FONT COLOR="#B22222">// * Real * Real = Real </FONT></I> ) * JxW[qp]; <I><FONT COLOR="#B22222">// ) * Real = Real </FONT></I> Me(i,j) += ( <I><FONT COLOR="#B22222">// ( </FONT></I> (1. - (dphase[qp] * dphase[qp])) <I><FONT COLOR="#B22222">// (Real - (Point * Point)) = Real </FONT></I> * phi[i][qp] * phi[j][qp] * weight[qp] <I><FONT COLOR="#B22222">// * Real * Real * Real = Real </FONT></I> ) * JxW[qp]; <I><FONT COLOR="#B22222">// ) * Real = Real </FONT></I> } <I><FONT COLOR="#B22222">// end of the matrix summation loop</FONT></I> } <I><FONT COLOR="#B22222">// end of quadrature point loop</FONT></I> Ke.add(1./speed , Ce); Ke.add(1./(speed*speed), Me); system.matrix->add_matrix (Ke, dof_indices); } <I><FONT COLOR="#B22222">// end of element loop</FONT></I> { <B><FONT COLOR="#228B22">const</FONT></B> <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> n_nodes = mesh.n_nodes(); <B><FONT COLOR="#A020F0">for</FONT></B> (<B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> n_cnt=0; n_cnt<n_nodes; n_cnt++) { <B><FONT COLOR="#228B22">const</FONT></B> Node& curr_node = mesh.node(n_cnt); <B><FONT COLOR="#A020F0">if</FONT></B> (fabs(curr_node(0)) < TOLERANCE && fabs(curr_node(1)) < TOLERANCE && fabs(curr_node(2)) < TOLERANCE) { <B><FONT COLOR="#228B22">unsigned</FONT></B> <B><FONT COLOR="#228B22">int</FONT></B> dn = curr_node.dof_number(0,0,0); system.rhs->add (dn, 1.); } } } #<B><FONT COLOR="#A020F0">else</FONT></B> assert(es.get_mesh().mesh_dimension() != 1); #endif <I><FONT COLOR="#B22222">//ifdef ENABLE_INFINITE_ELEMENTS</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 ./ex6-devel*************************************************************** ERROR: This example requires the library to be compiled with Infinite Element support![0] ex6.C, line 91, compiled Jun 6 2007 at 11:54:44 **************************************************************** Done Running Example ./ex6-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 + -