📄 ex14.php
字号:
<?php $root=""; ?><?php require($root."navigation.php"); ?><html><head> <?php load_style($root); ?></head> <body> <?php make_navigation("ex14",$root)?> <div class="content"><a name="comments"></a> <div class = "comment"><h1>Example 14 - Laplace Equation in the L-Shaped Domain</h1><br><br>This example solves the Laplace equation on the classic "L-shaped"domain with adaptive mesh refinement. In this case, the exactsolution is u(r,\theta) = r^{2/3} * \sin ( (2/3) * \theta), butthe standard Kelly error indicator is used to estimate the error.The initial mesh contains three QUAD9 elements which represent thestandard quadrants I, II, and III of the domain [-1,1]x[-1,1],i.e.Element 0: [-1,0]x[ 0,1]Element 1: [ 0,1]x[ 0,1]Element 2: [-1,0]x[-1,0]The mesh is provided in the standard libMesh ASCII format filenamed "lshaped.xda". In addition, an input file named "ex14.in"is provided which allows the user to set several parameters forthe solution so that the problem can be re-run without are-compile. The solution technique employed is to have arefinement loop with a linear solve inside followed by arefinement of the grid and projection of the solution to the new gridIn the final loop iteration, there is no additionalrefinement after the solve. In the input file "ex14.in", the variable"max_r_steps" controls the number of refinement steps,"max_r_level" controls the maximum element refinement level, and"refine_percentage" / "coarsen_percentage" determine the number ofelements which will be refined / coarsened at each step.<br><br>LibMesh include files.</div><div class ="fragment"><pre> #include "mesh.h" #include "equation_systems.h" #include "linear_implicit_system.h" #include "gmv_io.h" #include "tecplot_io.h" #include "fe.h" #include "quadrature_gauss.h" #include "dense_matrix.h" #include "dense_vector.h" #include "sparse_matrix.h" #include "mesh_refinement.h" #include "error_vector.h" #include "exact_error_estimator.h" #include "kelly_error_estimator.h" #include "patch_recovery_error_estimator.h" #include "uniform_refinement_estimator.h" #include "hp_coarsentest.h" #include "hp_singular.h" #include "mesh_generation.h" #include "mesh_modification.h" #include "getpot.h" #include "exact_solution.h" #include "dof_map.h" #include "numeric_vector.h" #include "elem.h" #include "string_to_enum.h" </pre></div><div class = "comment">Function prototype. This is the function that will assemblethe linear system for our Laplace problem. Note that thefunction will take the \p EquationSystems object and thename of the system we are assembling as input. From the\p EquationSystems object we have acess to the \p Mesh andother objects we might need.</div><div class ="fragment"><pre> void assemble_laplace(EquationSystems& es, const std::string& system_name); </pre></div><div class = "comment">Prototype for calculation of the exact solution. Usefulfor setting boundary conditions.</div><div class ="fragment"><pre> Number exact_solution(const Point& p, const Parameters&, // EquationSystem parameters, not needed const std::string&, // sys_name, not needed const std::string&); // unk_name, not needed); </pre></div><div class = "comment">Prototype for calculation of the gradient of the exact solution. </div><div class ="fragment"><pre> Gradient exact_derivative(const Point& p, const Parameters&, // EquationSystems parameters, not needed const std::string&, // sys_name, not needed const std::string&); // unk_name, not needed); </pre></div><div class = "comment">These are non-const because the input file may change it,It is global because our exact_* functions use it.<br><br>Set the dimensionality of the mesh</div><div class ="fragment"><pre> unsigned int dim = 2; </pre></div><div class = "comment">Choose whether or not to use the singular solution</div><div class ="fragment"><pre> bool singularity = true; int main(int argc, char** argv) {</pre></div><div class = "comment">Initialize libMesh.</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">Parse the input file</div><div class ="fragment"><pre> GetPot input_file("ex14.in"); </pre></div><div class = "comment">Read in parameters from the input file</div><div class ="fragment"><pre> const unsigned int max_r_steps = input_file("max_r_steps", 3); const unsigned int max_r_level = input_file("max_r_level", 3); const Real refine_percentage = input_file("refine_percentage", 0.5); const Real coarsen_percentage = input_file("coarsen_percentage", 0.5); const unsigned int uniform_refine = input_file("uniform_refine",0); const std::string refine_type = input_file("refinement_type", "h"); const std::string approx_type = input_file("approx_type", "LAGRANGE"); const unsigned int approx_order = input_file("approx_order", 1); const std::string element_type = input_file("element_type", "tensor"); const int extra_error_quadrature = input_file("extra_error_quadrature", 0); const int max_linear_iterations = input_file("max_linear_iterations", 5000); dim = input_file("dimension", 2); const std::string indicator_type = input_file("indicator_type", "kelly"); singularity = input_file("singularity", true); </pre></div><div class = "comment">Output file for plotting the error as a function ofthe number of degrees of freedom.</div><div class ="fragment"><pre> std::string approx_name = ""; if (element_type == "tensor") approx_name += "bi"; if (approx_order == 1) approx_name += "linear"; else if (approx_order == 2) approx_name += "quadratic"; else if (approx_order == 3) approx_name += "cubic"; else if (approx_order == 4) approx_name += "quartic"; std::string output_file = approx_name; output_file += "_"; output_file += refine_type; if (uniform_refine == 0) output_file += "_adaptive.m"; else output_file += "_uniform.m"; std::ofstream out (output_file.c_str()); out << "% dofs L2-error H1-error" << std::endl; out << "e = [" << std::endl; </pre></div><div class = "comment">Create an n-dimensional mesh.</div><div class ="fragment"><pre> Mesh mesh (dim); </pre></div><div class = "comment">Read in the mesh</div><div class ="fragment"><pre> if (dim == 1) MeshTools::Generation::build_line(mesh,1,-1.,0.); else if (dim == 2) mesh.read("lshaped.xda"); else mesh.read("lshaped3D.xda"); </pre></div><div class = "comment">Use triangles if the config file says so</div><div class ="fragment"><pre> if (element_type == "simplex") MeshTools::Modification::all_tri(mesh); </pre></div><div class = "comment">We used first order elements to describe the geometry,but we may need second order elements to hold the degreesof freedom</div><div class ="fragment"><pre> if (approx_order > 1 || refine_type != "h") mesh.all_second_order(); </pre></div><div class = "comment">Mesh Refinement object</div><div class ="fragment"><pre> MeshRefinement mesh_refinement(mesh); mesh_refinement.refine_fraction() = refine_percentage; mesh_refinement.coarsen_fraction() = coarsen_percentage; mesh_refinement.max_h_level() = max_r_level; </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 system and its variables.</div><div class ="fragment"><pre> {</pre></div><div class = "comment">Creates a system named "Laplace"</div><div class ="fragment"><pre> LinearImplicitSystem& system = equation_systems.add_system<LinearImplicitSystem> ("Laplace"); </pre></div><div class = "comment">Adds the variable "u" to "Laplace", using the finite element type and order specifiedin the config file</div><div class ="fragment"><pre> system.add_variable("u", static_cast<Order>(approx_order), Utility::string_to_enum<FEFamily>(approx_type)); </pre></div><div class = "comment">Give the system a pointer to the matrix assemblyfunction.</div><div class ="fragment"><pre> system.attach_assemble_function (assemble_laplace); </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">Set linear solver max iterations</div><div class ="fragment"><pre> equation_systems.parameters.set<unsigned int>("linear solver maximum iterations") = max_linear_iterations; </pre></div><div class = "comment">Linear solver tolerance.</div><div class ="fragment"><pre> equation_systems.parameters.set<Real>("linear solver tolerance") = TOLERANCE * TOLERANCE * TOLERANCE; </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">Construct ExactSolution object and attach solution functions</div><div class ="fragment"><pre> ExactSolution exact_sol(equation_systems); exact_sol.attach_exact_value(exact_solution);⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -