📄 ex0.php
字号:
<?php $root=""; ?><?php require($root."navigation.php"); ?><html><head> <?php load_style($root); ?></head> <body> <?php make_navigation("ex0",$root)?> <div class="content"><a name="comments"></a> <div class = "comment"><h1>Example 0 - Solving 1D PDE Using Adaptive Mesh Refinement</h1><br><br>This example demonstrates how to solve a simple 1D problemusing adaptive mesh refinement. The PDE that is solved is:-epsilon*u''(x) + u(x) = 1, on the domain [0,1] with boundary conditions u(0) = u(1) = 0 and where epsilon << 1.<br><br>The approach used to solve 1D problems in libMesh is virtually identical tosolving 2D or 3D problems, so in this sense this example represents a goodstarting point for new users. Note that many concepts are used in this example which are explained more fully in subsequent examples.<br><br>Libmesh includes</div><div class ="fragment"><pre> #include "mesh.h" #include "mesh_generation.h" #include "edge_edge3.h" #include "gnuplot_io.h" #include "equation_systems.h" #include "linear_implicit_system.h" #include "fe.h" #include "quadrature_gauss.h" #include "sparse_matrix.h" #include "dof_map.h" #include "numeric_vector.h" #include "dense_matrix.h" #include "dense_vector.h" #include "error_vector.h" #include "kelly_error_estimator.h" #include "mesh_refinement.h" void assemble_1D(EquationSystems& es, const std::string& system_name); int main(int argc, char** argv) { </pre></div><div class = "comment">Initialize the library. This is necessary because the librarymay depend on a number of other libraries (i.e. MPI and Petsc)that require initialization before use. </div><div class ="fragment"><pre> libMesh::init(argc, argv); #ifndef ENABLE_AMR std::cerr << "ERROR: This example requires libMesh to be\n" << "compiled with AMR support!" << std::endl; return 0; #else {</pre></div><div class = "comment">Create a new 1 dimensional mesh</div><div class ="fragment"><pre> const unsigned int dim = 1; Mesh mesh(dim); </pre></div><div class = "comment">Build a 1D mesh with 4 elements from x=0 to x=1, using EDGE3 (i.e. quadratic) 1D elements. They are called EDGE3 elementsbecause a quadratic element contains 3 nodes.</div><div class ="fragment"><pre> MeshTools::Generation::build_line(mesh,4,0.,1.,EDGE3); </pre></div><div class = "comment">Define the equation systems object and the system we are goingto solve. See Example 2 for more details.</div><div class ="fragment"><pre> EquationSystems equation_systems(mesh); LinearImplicitSystem& system = equation_systems.add_system <LinearImplicitSystem>("1D"); </pre></div><div class = "comment">Add a variable "u" to the system, using second-order approximation</div><div class ="fragment"><pre> system.add_variable("u",SECOND); </pre></div><div class = "comment">Give the system a pointer to the matrix assembly function. This will be called when needed by the library.</div><div class ="fragment"><pre> system.attach_assemble_function(assemble_1D); </pre></div><div class = "comment">Define the mesh refinement object that takes care of adaptivelyrefining the mesh.</div><div class ="fragment"><pre> MeshRefinement mesh_refinement(mesh); </pre></div><div class = "comment">These parameters determine the proportion of elements that willbe refined and coarsened. Any element within 30% of the maximum error on any element will be refined, and any element within 30% of the minimum error on any element might be coarsened</div><div class ="fragment"><pre> mesh_refinement.refine_fraction() = 0.7; mesh_refinement.coarsen_fraction() = 0.3;</pre></div><div class = "comment">We won't refine any element more than 5 times in total</div><div class ="fragment"><pre> mesh_refinement.max_h_level() = 5; </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">Refinement parameters</div><div class ="fragment"><pre> const unsigned int max_r_steps = 5; // Refine the mesh 5 times </pre></div><div class = "comment">Define the refinement loop</div><div class ="fragment"><pre> for(unsigned int r_step=0; r_step<=max_r_steps; r_step++) {</pre></div><div class = "comment">Solve the equation system</div><div class ="fragment"><pre> equation_systems.get_system("1D").solve(); </pre></div><div class = "comment">We need to ensure that the mesh is not refined on the last iterationof this loop, since we do not want to refine the mesh unless we aregoing to solve the equation system for that refined mesh.</div><div class ="fragment"><pre> if(r_step != max_r_steps) {</pre></div><div class = "comment">Define object for error estimation, see Example 10 for more details.</div><div class ="fragment"><pre> ErrorVector error; KellyErrorEstimator error_estimator; </pre></div><div class = "comment">Compute the error for each active element</div><div class ="fragment"><pre> error_estimator.estimate_error(system, error); </pre></div><div class = "comment">Flag elements to be refined and coarsened</div><div class ="fragment"><pre> mesh_refinement.flag_elements_by_error_fraction (error); </pre></div><div class = "comment">Perform refinement and coarsening</div><div class ="fragment"><pre> mesh_refinement.refine_and_coarsen_elements(); </pre></div><div class = "comment">Reinitialize the equation_systems object for the newly refinedmesh. One of the steps in this is project the solution onto the new mesh</div><div class ="fragment"><pre> equation_systems.reinit(); } } </pre></div><div class = "comment">Construct gnuplot plotting object, pass in mesh, title of plotand boolean to indicate use of grid in plot. The grid is used toshow the edges of each element in the mesh.</div><div class ="fragment"><pre> GnuPlotIO plot(mesh,"Example 0", GnuPlotIO::GRID_ON); </pre></div><div class = "comment">Write out script to be called from within gnuplot:Load gnuplot, then type "call 'gnuplot_script'" from gnuplot prompt</div><div class ="fragment"><pre> plot.write_equation_systems("gnuplot_script",equation_systems); } #endif // #ifndef ENABLE_AMR </pre></div><div class = "comment">All done. Call the libMesh::close() function to close anyexternal libraries and check for leaked memory. To be absoluteycertain this is called last we will return its value. Thisalso allows main to return nonzero if memory is leaked, whichcan be useful for testing purposes.</div><div class ="fragment"><pre> return libMesh::close(); } </pre></div><div class = "comment">Define the matrix assembly function for the 1D PDE we are solving</div><div class ="fragment"><pre> void assemble_1D(EquationSystems& es, const std::string& system_name) { #ifdef ENABLE_AMR </pre></div><div class = "comment">It is a good idea to check we are solving the correct system</div><div class ="fragment"><pre> assert(system_name == "1D"); </pre></div><div class = "comment">Get a reference to the mesh object</div><div class ="fragment"><pre> const Mesh& mesh = es.get_mesh(); </pre></div><div class = "comment">The dimension we are using, i.e. dim==1</div><div class ="fragment"><pre> const unsigned int dim = mesh.mesh_dimension(); </pre></div><div class = "comment">Get a reference to the system we are solving</div><div class ="fragment"><pre> LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("1D"); </pre></div><div class = "comment">Get a reference to the DofMap object for this system. The DofMap objecthandles the index translation from node and element numbers to degree offreedom numbers. DofMap's are discussed in more detail in future examples.</div><div class ="fragment"><pre> const DofMap& dof_map = system.get_dof_map(); </pre></div><div class = "comment">Get a constant reference to the Finite Element type for the first (and only) variable in the system.</div><div class ="fragment"><pre> FEType fe_type = dof_map.variable_type(0); </pre></div><div class = "comment">Build a finite element object of the specified type. The buildfunction dynamically allocates memory so we use an AutoPtr in this case.An AutoPtr is a pointer that cleans up after itself. See examples 3 and 4for more details on AutoPtr.</div><div class ="fragment"><pre> AutoPtr<FEBase> fe(FEBase::build(dim, fe_type)); </pre></div><div class = "comment">Tell the finite element object to use fifth order Gaussian quadrature</div><div class ="fragment"><pre> QGauss qrule(dim,FIFTH); fe->attach_quadrature_rule(&qrule); </pre></div><div class = "comment">Here we define some references to cell-specific data that will be used to assemble the linear system.<br><br>The element Jacobian * quadrature weight at each integration point.</div><div class ="fragment"><pre> const std::vector<Real>& JxW = fe->get_JxW(); </pre></div>⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -