📄 ex17.php
字号:
<?php $root=""; ?><?php require($root."navigation.php"); ?><html><head> <?php load_style($root); ?></head> <body> <?php make_navigation("ex17",$root)?> <div class="content"><a name="comments"></a> <div class = "comment"></div><div class ="fragment"><pre> #include "libmesh.h" #include "mesh.h" #include "mesh_generation.h" #include "gmv_io.h" #include "eigen_system.h" #include "equation_systems.h" #include "fe.h" #include "quadrature_gauss.h" #include "dense_matrix.h" #include "sparse_matrix.h" #include "numeric_vector.h" #include "dof_map.h" </pre></div><div class = "comment">Function prototype. This is the function that will assemblethe eigen system. Here, we will simply assemble a mass matrix.</div><div class ="fragment"><pre> void assemble_mass(EquationSystems& es, const std::string& system_name); int main (int argc, char** argv) {</pre></div><div class = "comment">Initialize libMesh and the dependent libraries.</div><div class ="fragment"><pre> libMesh::init (argc, argv); </pre></div><div class = "comment">This example is designed for the SLEPc eigen solver interface.</div><div class ="fragment"><pre> #ifndef HAVE_SLEPC std::cerr << "ERROR: This example requires libMesh to be\n" << "compiled with SLEPc eigen solvers support!" << std::endl; return 0; #else </pre></div><div class = "comment">Braces are used to force object scope. </div><div class ="fragment"><pre> {</pre></div><div class = "comment">Check for proper usage.</div><div class ="fragment"><pre> if (argc < 3) { std::cerr << "\nUsage: " << argv[0] << " -n <number of eigen values>" << std::endl; error(); } </pre></div><div class = "comment">Tell the user what we are doing.</div><div class ="fragment"><pre> else { std::cout << "Running " << argv[0]; for (int i=1; i<argc; i++) std::cout << " " << argv[i]; std::cout << std::endl << std::endl; } </pre></div><div class = "comment">Set the dimensionality.</div><div class ="fragment"><pre> const unsigned int dim = 2; </pre></div><div class = "comment">Get the number of eigen values to be computed from argv[2]</div><div class ="fragment"><pre> const unsigned int nev = std::atoi(argv[2]); </pre></div><div class = "comment">Create a dim-dimensional mesh.</div><div class ="fragment"><pre> Mesh mesh (dim); </pre></div><div class = "comment">Use the internal mesh generator to create a uniformgrid on a square.</div><div class ="fragment"><pre> MeshTools::Generation::build_square (mesh, 20, 20, -1., 1., -1., 1., QUAD4); </pre></div><div class = "comment">Print information about the mesh to the screen.</div><div class ="fragment"><pre> mesh.print_info(); </pre></div><div class = "comment">Create an equation systems object.</div><div class ="fragment"><pre> EquationSystems equation_systems (mesh); </pre></div><div class = "comment">Create a EigenSystem named "Eigensystem" and (for convenience)use a reference to the system we create.</div><div class ="fragment"><pre> EigenSystem & eigen_system = equation_systems.add_system<EigenSystem> ("Eigensystem"); </pre></div><div class = "comment">Declare the system variables.</div><div class ="fragment"><pre> {</pre></div><div class = "comment">Adds the variable "p" to "Eigensystem". "p"will be approximated using second-order approximation.</div><div class ="fragment"><pre> eigen_system.add_variable("p", FIRST); </pre></div><div class = "comment">Give the system a pointer to the matrix assemblyfunction defined below.</div><div class ="fragment"><pre> eigen_system.attach_assemble_function (assemble_mass); </pre></div><div class = "comment">Set necessary parametrs used in EigenSystem::solve(),i.e. the number of requested eigenpairs \p nev and the numberof basis vectors \p ncv used in the solution algorithm. Note thatncv >= nev must hold and ncv >= 2*nev is recommended.</div><div class ="fragment"><pre> equation_systems.parameters.set<unsigned int>("eigenpairs") = nev; equation_systems.parameters.set<unsigned int>("basis vectors") = nev*3; </pre></div><div class = "comment">Set the eigen solver type. SLEPc offers various solvers such asthe Arnoldi and subspace method. Italso offers interfaces to other solver packages (e.g. ARPACK).</div><div class ="fragment"><pre> eigen_system.eigen_solver->set_eigensolver_type(ARNOLDI); </pre></div><div class = "comment">Set the solver tolerance and the maximum number of iterations. </div><div class ="fragment"><pre> equation_systems.parameters.set<Real>("linear solver tolerance") = pow(TOLERANCE, 5./3.); equation_systems.parameters.set<unsigned int> ("linear solver maximum iterations") = 1000; </pre></div><div class = "comment">Set the type of the problem, here we deal witha generalized Hermitian problem.</div><div class ="fragment"><pre> eigen_system.set_eigenproblem_type(GHEP); </pre></div><div class = "comment">Set the eigenvalues to be computed. Note that notall solvers support this.eigen_system.eigen_solver->set_position_of_spectrum(SMALLEST_MAGNITUDE);<br><br>Initialize the data structures for the equation system.</div><div class ="fragment"><pre> equation_systems.init(); </pre></div><div class = "comment">Prints information about the system to the screen.</div><div class ="fragment"><pre> equation_systems.print_info(); } </pre></div><div class = "comment">Solve the system "Eigensystem".</div><div class ="fragment"><pre> eigen_system.solve(); </pre></div><div class = "comment">Get the number of converged eigen pairs.</div><div class ="fragment"><pre> unsigned int nconv = eigen_system.get_n_converged(); std::cout << "Number of converged eigenpairs: " << nconv << "\n" << std::endl; </pre></div><div class = "comment">Get the last converged eigenpair</div><div class ="fragment"><pre> if (nconv != 0) { eigen_system.get_eigenpair(nconv-1); </pre></div><div class = "comment">Write the eigen vector to file.</div><div class ="fragment"><pre> char buf[14]; sprintf (buf, "out.gmv"); GMVIO (mesh).write_equation_systems (buf, equation_systems); } else { std::cout << "WARNING: Solver did not converge!\n" << nconv << std::endl; } } #endif // HAVE_SLEPC </pre></div><div class = "comment">All done. </div><div class ="fragment"><pre> return libMesh::close (); } void assemble_mass(EquationSystems& es, const std::string& system_name) { </pre></div><div class = "comment">It is a good idea to make sure we are assemblingthe proper system.</div><div class ="fragment"><pre> assert (system_name == "Eigensystem"); #ifdef HAVE_SLEPC </pre></div><div class = "comment">Get a constant reference to the mesh object.</div><div class ="fragment"><pre> const Mesh& mesh = es.get_mesh(); </pre></div><div class = "comment">The dimension that we are running.</div><div class ="fragment"><pre> const unsigned int dim = mesh.mesh_dimension(); </pre></div><div class = "comment">Get a reference to our system.</div><div class ="fragment"><pre> EigenSystem & eigen_system = es.get_system<EigenSystem> (system_name); </pre></div><div class = "comment">Get a constant reference to the Finite Element typefor the first (and only) variable in the system.</div><div class ="fragment"><pre> FEType fe_type = eigen_system.get_dof_map().variable_type(0); </pre></div><div class = "comment">A reference to the two system matrices</div><div class ="fragment"><pre> SparseMatrix<Number>& matrix_A = *eigen_system.matrix_A; SparseMatrix<Number>& matrix_B = *eigen_system.matrix_B; </pre></div><div class = "comment">Build a Finite Element object of the specified type. Since the\p FEBase::build() member dynamically creates memory we willstore the object as an \p AutoPtr<FEBase>. This can be thoughtof as a pointer that will clean up after itself.</div><div class ="fragment"><pre> AutoPtr<FEBase> fe (FEBase::build(dim, fe_type)); </pre></div><div class = "comment">A Gauss quadrature rule for numerical integration.Use the default quadrature order.</div><div class ="fragment"><pre> QGauss qrule (dim, fe_type.default_quadrature_order()); </pre></div><div class = "comment">Tell the finite element object to use our quadrature rule.</div><div class ="fragment"><pre> fe->attach_quadrature_rule (&qrule); </pre></div><div class = "comment">The element Jacobian * quadrature weight at each integration point. ⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -