ex14.php

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

            exact_sol.attach_exact_deriv(exact_derivative);
        
</pre>
</div>
<div class = "comment">
Use higher quadrature order for more accurate error results
</div>

<div class ="fragment">
<pre>
            exact_sol.extra_quadrature_order(extra_error_quadrature);
        
</pre>
</div>
<div class = "comment">
Convenient reference to the system
</div>

<div class ="fragment">
<pre>
            LinearImplicitSystem& system =
              equation_systems.get_system&lt;LinearImplicitSystem&gt;("Laplace");
        
</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++)
              {
        	std::cout &lt;&lt; "Beginning Solve " &lt;&lt; r_step &lt;&lt; std::endl;
        	
</pre>
</div>
<div class = "comment">
Solve the system "Laplace", just like example 2.
</div>

<div class ="fragment">
<pre>
                system.solve();
        
        	std::cout &lt;&lt; "System has: " &lt;&lt; equation_systems.n_active_dofs()
        		  &lt;&lt; " degrees of freedom."
        		  &lt;&lt; std::endl;
        
        	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("Laplace", "u");
        
</pre>
</div>
<div class = "comment">
Print out the error values
</div>

<div class ="fragment">
<pre>
                std::cout &lt;&lt; "L2-Error is: "
        		  &lt;&lt; exact_sol.l2_error("Laplace", "u")
        		  &lt;&lt; std::endl;
        	std::cout &lt;&lt; "H1-Error is: "
        		  &lt;&lt; exact_sol.h1_error("Laplace", "u")
        		  &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("Laplace", "u") &lt;&lt; " "
        	    &lt;&lt; exact_sol.h1_error("Laplace", "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)
                      {
        
</pre>
</div>
<div class = "comment">
The \p ErrorVector is a particular \p StatisticsVector
for computing error information on a finite element mesh.
</div>

<div class ="fragment">
<pre>
                        ErrorVector error;
        		
                        if (indicator_type == "exact")
                          {
</pre>
</div>
<div class = "comment">
The \p ErrorEstimator class interrogates a
finite element solution and assigns to each
element a positive error value.
This value is used for deciding which elements to
refine and which to coarsen.
For these simple test problems, we can use
numerical quadrature of the exact error between
the approximate and analytic solutions.
However, for real problems, we would need an error
indicator which only relies on the approximate
solution.
</div>

<div class ="fragment">
<pre>
                            ExactErrorEstimator error_estimator;
        
                            error_estimator.attach_exact_value(exact_solution);
                            error_estimator.attach_exact_deriv(exact_derivative);
        
</pre>
</div>
<div class = "comment">
We optimize in H1 norm
</div>

<div class ="fragment">
<pre>
                            error_estimator.sobolev_order() = 1;
        
</pre>
</div>
<div class = "comment">
Compute the error for each active element using
the provided indicator.  Note in general you
will need to provide an error estimator
specifically designed for your application.
</div>

<div class ="fragment">
<pre>
                            error_estimator.estimate_error (system, error);
                          }
                        else if (indicator_type == "patch")
                          {
</pre>
</div>
<div class = "comment">
The patch recovery estimator should give a
good estimate of the solution interpolation
error.
</div>

<div class ="fragment">
<pre>
                            PatchRecoveryErrorEstimator error_estimator;
        
        		    error_estimator.estimate_error (system, error);
                          }
                        else if (indicator_type == "uniform")
                          {
</pre>
</div>
<div class = "comment">
Error indication based on uniform refinement
is reliable, but very expensive.
</div>

<div class ="fragment">
<pre>
                            UniformRefinementEstimator error_estimator;
        
        		    error_estimator.estimate_error (system, error);
                          }
                        else
                          {
                            assert (indicator_type == "kelly");
        
</pre>
</div>
<div class = "comment">
The Kelly error estimator is based on 
an error bound for the Poisson problem
on linear elements, but is useful for
driving adaptive refinement in many problems
</div>

<div class ="fragment">
<pre>
                            KellyErrorEstimator error_estimator;
        
        		    error_estimator.estimate_error (system, error);
                          }
        		
</pre>
</div>
<div class = "comment">
This takes the error in \p error and decides which elements
will be coarsened or refined.  Any element within 20% of the
maximum error on any element will be refined, and any
element within 10% of the minimum error on any element might
be coarsened. Note that the elements flagged for refinement
will be refined, but those flagged for coarsening _might_ be
coarsened.
</div>

<div class ="fragment">
<pre>
                        mesh_refinement.flag_elements_by_error_fraction (error);
        
</pre>
</div>
<div class = "comment">
If we are doing adaptive p refinement, we want
elements flagged for that instead.
</div>

<div class ="fragment">
<pre>
                        if (refine_type == "p")
                          mesh_refinement.switch_h_to_p_refinement();
</pre>
</div>
<div class = "comment">
If we are doing "matched hp" refinement, we
flag elements for both h and p
</div>

<div class ="fragment">
<pre>
                        if (refine_type == "matchedhp")
                          mesh_refinement.add_p_to_h_refinement();
</pre>
</div>
<div class = "comment">
If we are doing hp refinement, we 
try switching some elements from h to p
</div>

<div class ="fragment">
<pre>
                        if (refine_type == "hp")
        	          {
        		    HPCoarsenTest hpselector;
                            hpselector.select_refinement(system);
        	          }
</pre>
</div>
<div class = "comment">
If we are doing "singular hp" refinement, we 
try switching most elements from h to p
</div>

<div class ="fragment">
<pre>
                        if (refine_type == "singularhp")
        	          {
</pre>
</div>
<div class = "comment">
This only differs from p refinement for
the singular problem
</div>

<div class ="fragment">
<pre>
                            assert (singularity);
        		    HPSingularity hpselector;
</pre>
</div>
<div class = "comment">
Our only singular point is at the origin
</div>

<div class ="fragment">
<pre>
                            hpselector.singular_points.push_back(Point());
                            hpselector.select_refinement(system);
        	          }
        		
</pre>
</div>
<div class = "comment">
This call actually refines and coarsens the flagged
elements.
</div>

<div class ="fragment">
<pre>
                        mesh_refinement.refine_and_coarsen_elements();
        	      }
        
        	    else if (uniform_refine == 1)
                      {
                        if (refine_type == "h" || refine_type == "hp" ||
                            refine_type == "matchedhp")
                          mesh_refinement.uniformly_refine(1);
                        if (refine_type == "p" || refine_type == "hp" ||
                            refine_type == "matchedhp")
                          mesh_refinement.uniformly_p_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 ("lshaped.gmv",
            					 equation_systems);
        
</pre>
</div>
<div class = "comment">
Close up the output file.
</div>

<div class ="fragment">
<pre>
            out &lt;&lt; "];" &lt;&lt; std::endl;
            out &lt;&lt; "hold on" &lt;&lt; std::endl;
            out &lt;&lt; "plot(e(:,1), e(:,2), 'bo-');" &lt;&lt; std::endl;
            out &lt;&lt; "plot(e(:,1), e(:,3), 'ro-');" &lt;&lt; std::endl;
</pre>
</div>
<div class = "comment">
out << "set(gca,'XScale', 'Log');" << std::endl;
out << "set(gca,'YScale', 'Log');" << std::endl;
</div>

<div class ="fragment">
<pre>
            out &lt;&lt; "xlabel('dofs');" &lt;&lt; std::endl;
            out &lt;&lt; "title('" &lt;&lt; approx_name &lt;&lt; " elements');" &lt;&lt; std::endl;
            out &lt;&lt; "legend('L2-error', 'H1-error');" &lt;&lt; std::endl;
</pre>
</div>
<div class = "comment">
out << "disp('L2-error linear fit');" << std::endl;
out << "polyfit(log10(e(:,1)), log10(e(:,2)), 1)" << std::endl;
out << "disp('H1-error linear fit');" << std::endl;
out << "polyfit(log10(e(:,1)), log10(e(:,3)), 1)" << std::endl;
</div>

<div class ="fragment">
<pre>
          }
        #endif // #ifndef ENABLE_AMR
          
</pre>
</div>
<div class = "comment">