📄 ex3.php
字号:
<?php $root=""; ?><?php require($root."navigation.php"); ?><html><head> <?php load_style($root); ?></head> <body> <?php make_navigation("ex3",$root)?> <div class="content"><a name="comments"></a> <div class = "comment"><h1>Example 3 - Solving a Poisson Problem</h1><br><br>This is the third example program. It builds onthe second example program by showing how to solve a simplePoisson system. This example also introduces the notionof customized matrix assembly functions, working with anexact solution, and using element iterators.We will not comment on things thatwere already explained in the second example.<br><br>C++ include files that we need</div><div class ="fragment"><pre> #include <iostream> #include <algorithm> #include <math.h> </pre></div><div class = "comment">Basic include files needed for the mesh functionality.</div><div class ="fragment"><pre> #include "libmesh.h" #include "mesh.h" #include "mesh_generation.h" #include "gmv_io.h" #include "linear_implicit_system.h" #include "equation_systems.h" </pre></div><div class = "comment">Define the Finite Element object.</div><div class ="fragment"><pre> #include "fe.h" </pre></div><div class = "comment">Define Gauss quadrature rules.</div><div class ="fragment"><pre> #include "quadrature_gauss.h" </pre></div><div class = "comment">Define useful datatypes for finite elementmatrix and vector components.</div><div class ="fragment"><pre> #include "sparse_matrix.h" #include "numeric_vector.h" #include "dense_matrix.h" #include "dense_vector.h" #include "elem.h" </pre></div><div class = "comment">Define the DofMap, which handles degree of freedomindexing.</div><div class ="fragment"><pre> #include "dof_map.h" #include "boundary_mesh.h" #include "boundary_info.h" </pre></div><div class = "comment">Function prototype. This is the function that will assemblethe linear system for our Poisson problem. Note that thefunction will take the EquationSystems object and thename of the system we are assembling as input. From theEquationSystems object we have access to the Mesh andother objects we might need.</div><div class ="fragment"><pre> void assemble_poisson(EquationSystems& es, const std::string& system_name); </pre></div><div class = "comment">Function prototype for the exact solution.</div><div class ="fragment"><pre> Real exact_solution (const Real x, const Real y, const Real z = 0.); int main (int argc, char** argv) { </pre></div><div class = "comment">Initialize Petsc, like in example 2.</div><div class ="fragment"><pre> libMesh::init (argc, argv); </pre></div><div class = "comment">Braces are used to force object scope, like in example 2</div><div class ="fragment"><pre> { </pre></div><div class = "comment">Brief message to the user regarding the program nameand command line arguments.</div><div class ="fragment"><pre> 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">Create a 2D mesh.</div><div class ="fragment"><pre> Mesh mesh (2); </pre></div><div class = "comment">Use the MeshTools::Generation mesh generator to create a uniformgrid on the square [-1,1]^2. We instruct the mesh generatorto build a mesh of 15x15 QUAD9 elements. Building QUAD9elements instead of the default QUAD4's we used in example 2allow us to use higher-order approximation.</div><div class ="fragment"><pre> MeshTools::Generation::build_square (mesh, 15, 15, -1., 1., -1., 1., QUAD9); </pre></div><div class = "comment">Print information about the mesh to the screen.Note that 5x5 QUAD9 elements actually has 11x11 nodes,so this mesh is significantly larger than the one in example 2.</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">Declare the Poisson system and its variables.The Poisson system is another example of a steady system.</div><div class ="fragment"><pre> equation_systems.add_system<LinearImplicitSystem> ("Poisson"); </pre></div><div class = "comment">Adds the variable "u" to "Poisson". "u"will be approximated using second-order approximation.</div><div class ="fragment"><pre> equation_systems.get_system("Poisson").add_variable("u", SECOND); </pre></div><div class = "comment">Give the system a pointer to the matrix assemblyfunction. This will be called when needed by thelibrary.</div><div class ="fragment"><pre> equation_systems.get_system("Poisson").attach_assemble_function (assemble_poisson); </pre></div><div class = "comment">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 "Poisson". Note that calling thismember will assemble the linear system and invokethe default Petsc solver, however the solver can becontrolled from the command line. For example,you can invoke conjugate gradient with:<br><br>./ex3 -ksp_type cg<br><br>and you can get a nice X-window that monitors the solverconvergence with:<br><br>./ex3 -ksp_xmonitor<br><br>if you linked against the appropriate X libraries when youbuilt PETSc.</div><div class ="fragment"><pre> equation_systems.get_system("Poisson").solve(); </pre></div><div class = "comment">After solving the system write the solutionto a GMV-formatted plot file.</div><div class ="fragment"><pre> GMVIO (mesh).write_equation_systems ("out.gmv", equation_systems); for (unsigned int i=0; i<3; i++) { here(); equation_systems.get_system("Poisson").solve(); BoundaryMesh boundary_mesh (mesh.mesh_dimension()-1); mesh.boundary_info->sync(boundary_mesh, false);</pre></div><div class = "comment">GMVIO(boundary_mesh).write("boundary.gmv");</div><div class ="fragment"><pre> } } </pre></div><div class = "comment">All done. </div><div class ="fragment"><pre> return libMesh::close(); } </pre></div><div class = "comment">We now define the matrix assembly function for thePoisson system. We need to first compute elementmatrices and right-hand sides, and then take intoaccount the boundary conditions, which will be handledvia a penalty method.</div><div class ="fragment"><pre> void assemble_poisson(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 == "Poisson"); </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 the LinearImplicitSystem we are solving</div><div class ="fragment"><pre> LinearImplicitSystem& system = es.get_system<LinearImplicitSystem> ("Poisson"); </pre></div><div class = "comment">A reference to the DofMap object for this system. The DofMapobject handles the index translation from node and element numbersto degree of freedom numbers. We will talk more about the DofMapin future examples.</div>⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -