ex15.php

一个用来实现偏微分方程中网格的计算库

<div class ="fragment">
<pre>
              equation_systems.parameters.set&lt;unsigned int&gt;
        		      ("linear solver maximum iterations") =
                              max_linear_iterations;
        
</pre>
</div>
<div class = "comment">
Linear solver tolerance.
</div>

<div class ="fragment">
<pre>
              equation_systems.parameters.set&lt;Real&gt;
        		      ("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 function to compute exact solution
</div>

<div class ="fragment">
<pre>
            ExactSolution exact_sol(equation_systems);
            exact_sol.attach_exact_value(exact_solution);
            exact_sol.attach_exact_deriv(exact_derivative);
            exact_sol.attach_exact_hessian(exact_hessian);
        
</pre>
</div>
<div class = "comment">
Construct zero solution object, useful for computing solution norms
Attaching "zero_solution" functions is unnecessary
</div>

<div class ="fragment">
<pre>
            ExactSolution zero_sol(equation_systems);
        
</pre>
</div>
<div class = "comment">
Convenient reference to the system
</div>

<div class ="fragment">
<pre>
            LinearImplicitSystem& system = 
              equation_systems.get_system&lt;LinearImplicitSystem&gt;("Biharmonic");
        
</pre>
</div>
<div class = "comment">
A refinement loop.
</div>

<div class ="fragment">
<pre>
            for (unsigned int r_step=0; r_step&lt;max_r_steps; r_step++)
              {
                mesh.print_info();
                equation_systems.print_info();
        
        	std::cout &lt;&lt; "Beginning Solve " &lt;&lt; r_step &lt;&lt; std::endl;
        	
</pre>
</div>
<div class = "comment">
Solve the system "Biharmonic", just like example 2.
</div>

<div class ="fragment">
<pre>
                system.solve();
        
        	std::cout &lt;&lt; "Linear solver converged at step: "
        		  &lt;&lt; system.n_linear_iterations()
        		  &lt;&lt; ", final residual: "
        		  &lt;&lt; system.final_linear_residual()
        		  &lt;&lt; std::endl;
        
</pre>
</div>
<div class = "comment">
Compute the error.
</div>

<div class ="fragment">
<pre>
                exact_sol.compute_error("Biharmonic", "u");
</pre>
</div>
<div class = "comment">
Compute the norm.
</div>

<div class ="fragment">
<pre>
                zero_sol.compute_error("Biharmonic", "u");
        
</pre>
</div>
<div class = "comment">
Print out the error values
</div>

<div class ="fragment">
<pre>
                std::cout &lt;&lt; "L2-Norm is: "
        		  &lt;&lt; zero_sol.l2_error("Biharmonic", "u")
        		  &lt;&lt; std::endl;
        	std::cout &lt;&lt; "H1-Norm is: "
        		  &lt;&lt; zero_sol.h1_error("Biharmonic", "u")
        		  &lt;&lt; std::endl;
        	std::cout &lt;&lt; "H2-Norm is: "
        		  &lt;&lt; zero_sol.h2_error("Biharmonic", "u")
        		  &lt;&lt; std::endl
        		  &lt;&lt; std::endl;
        	std::cout &lt;&lt; "L2-Error is: "
        		  &lt;&lt; exact_sol.l2_error("Biharmonic", "u")
        		  &lt;&lt; std::endl;
        	std::cout &lt;&lt; "H1-Error is: "
        		  &lt;&lt; exact_sol.h1_error("Biharmonic", "u")
        		  &lt;&lt; std::endl;
        	std::cout &lt;&lt; "H2-Error is: "
        		  &lt;&lt; exact_sol.h2_error("Biharmonic", "u")
        		  &lt;&lt; std::endl
        		  &lt;&lt; std::endl;
        
</pre>
</div>
<div class = "comment">
Print to output file
</div>

<div class ="fragment">
<pre>
                out &lt;&lt; equation_systems.n_active_dofs() &lt;&lt; " "
        	    &lt;&lt; exact_sol.l2_error("Biharmonic", "u") &lt;&lt; " "
        	    &lt;&lt; exact_sol.h1_error("Biharmonic", "u") &lt;&lt; " "
        	    &lt;&lt; exact_sol.h2_error("Biharmonic", "u") &lt;&lt; std::endl;
        
</pre>
</div>
<div class = "comment">
Possibly refine the mesh
</div>

<div class ="fragment">
<pre>
                if (r_step+1 != max_r_steps)
        	  {
        	    std::cout &lt;&lt; "  Refining the mesh..." &lt;&lt; std::endl;
        
        	    if (uniform_refine == 0)
        	      {
        		ErrorVector error;
        		LaplacianErrorEstimator error_estimator;
        
        		error_estimator.estimate_error(system, error);
                        mesh_refinement.flag_elements_by_elem_fraction (error);
        
        		std::cerr &lt;&lt; "Mean Error: " &lt;&lt; error.mean() &lt;&lt;
        				std::endl;
        		std::cerr &lt;&lt; "Error Variance: " &lt;&lt; error.variance() &lt;&lt;
        				std::endl;
        
        		mesh_refinement.refine_and_coarsen_elements();
                      }
        	    else
        	      {
                        mesh_refinement.uniformly_refine(1);
                      }
        		
</pre>
</div>
<div class = "comment">
This call reinitializes the \p EquationSystems object for
the newly refined mesh.  One of the steps in the
reinitialization is projecting the \p solution,
\p old_solution, etc... vectors from the old mesh to
the current one.
</div>

<div class ="fragment">
<pre>
                    equation_systems.reinit ();
        	  }
              }	    
            
</pre>
</div>
<div class = "comment">
Write out the solution
After solving the system write the solution
to a GMV-formatted plot file.
</div>

<div class ="fragment">
<pre>
            GMVIO (mesh).write_equation_systems (gmv_file,
            					 equation_systems);
</pre>
</div>
<div class = "comment">
Close up the output file.
</div>

<div class ="fragment">
<pre>
            out &lt;&lt; "];\n"
        	&lt;&lt; "hold on\n"
        	&lt;&lt; "loglog(e(:,1), e(:,2), 'bo-');\n"
        	&lt;&lt; "loglog(e(:,1), e(:,3), 'ro-');\n"
        	&lt;&lt; "loglog(e(:,1), e(:,4), 'go-');\n"
        	&lt;&lt; "xlabel('log(dofs)');\n"
        	&lt;&lt; "ylabel('log(error)');\n"
        	&lt;&lt; "title('C1 " &lt;&lt; approx_type &lt;&lt; " elements');\n"
        	&lt;&lt; "legend('L2-error', 'H1-error', 'H2-error');\n";
          }
          
</pre>
</div>
<div class = "comment">
All done.  
</div>

<div class ="fragment">
<pre>
          return libMesh::close ();
        #endif // #ifndef ENABLE_SECOND_DERIVATIVES
        #endif // #ifndef ENABLE_AMR
        }
        
        
        
        Number exact_2D_solution(const Point& p,
        		         const Parameters&,  // parameters, not needed
        		         const std::string&, // sys_name, not needed
        		         const std::string&) // unk_name, not needed
        {
</pre>
</div>
<div class = "comment">
x and y coordinates in space
</div>

<div class ="fragment">
<pre>
          const Real x = p(0);
          const Real y = p(1);
        
</pre>
</div>
<div class = "comment">
analytic solution value
</div>

<div class ="fragment">
<pre>
          return 256.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y);
        }
        
        
</pre>
</div>
<div class = "comment">
We now define the gradient of the exact solution
</div>

<div class ="fragment">
<pre>
        Gradient exact_2D_derivative(const Point& p,
        			     const Parameters&,  // parameters, not needed
        			     const std::string&, // sys_name, not needed
        			     const std::string&) // unk_name, not needed
        {
</pre>
</div>
<div class = "comment">
x and y coordinates in space
</div>

<div class ="fragment">
<pre>
          const Real x = p(0);
          const Real y = p(1);
        
</pre>
</div>
<div class = "comment">
First derivatives to be returned.
</div>

<div class ="fragment">
<pre>
          Gradient gradu;
        
          gradu(0) = 256.*2.*(x-x*x)*(1-2*x)*(y-y*y)*(y-y*y);
          gradu(1) = 256.*2.*(x-x*x)*(x-x*x)*(y-y*y)*(1-2*y);
        
          return gradu;
        }
        
        
</pre>
</div>
<div class = "comment">
We now define the hessian of the exact solution
</div>

<div class ="fragment">
<pre>
        Tensor exact_2D_hessian(const Point& p,
        			const Parameters&,  // parameters, not needed
        			const std::string&, // sys_name, not needed
        			const std::string&) // unk_name, not needed
        {
</pre>
</div>
<div class = "comment">
Second derivatives to be returned.
</div>

<div class ="fragment">
<pre>
          Tensor hessu;
          
</pre>
</div>
<div class = "comment">
x and y coordinates in space
</div>

<div class ="fragment">
<pre>
          const Real x = p(0);
          const Real y = p(1);
        
          hessu(0,0) = 256.*2.*(1-6.*x+6.*x*x)*(y-y*y)*(y-y*y);
          hessu(0,1) = 256.*4.*(x-x*x)*(1.-2.*x)*(y-y*y)*(1.-2.*y);
          hessu(1,1) = 256.*2.*(x-x*x)*(x-x*x)*(1.-6.*y+6.*y*y);
        
</pre>
</div>
<div class = "comment">
Hessians are always symmetric
</div>

<div class ="fragment">
<pre>
          hessu(1,0) = hessu(0,1);
          return hessu;
        }
        
        
        
        Number forcing_function_2D(const Point& p)
        {
</pre>
</div>
<div class = "comment">
x and y coordinates in space
</div>

<div class ="fragment">
<pre>
          const Real x = p(0);
          const Real y = p(1);
        
</pre>
</div>
<div class = "comment">
Equals laplacian(laplacian(u))
</div>

<div class ="fragment">
<pre>
          return 256. * 8. * (3.*((y-y*y)*(y-y*y)+(x-x*x)*(x-x*x))
                 + (1.-6.*x+6.*x*x)*(1.-6.*y+6.*y*y));
        }
        
        
        
        Number exact_3D_solution(const Point& p,
        		         const Parameters&,  // parameters, not needed
        		         const std::string&, // sys_name, not needed
        		         const std::string&) // unk_name, not needed
        {