ex14.c

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

/* $Id: ex14.C 2837 2008-05-08 17:23:37Z roystgnr $ */

/* The Next Great Finite Element Library. */
/* Copyright (C) 2004  Benjamin S. Kirk, John W. Peterson */

/* This library is free software; you can redistribute it and/or */
/* modify it under the terms of the GNU Lesser General Public */
/* License as published by the Free Software Foundation; either */
/* version 2.1 of the License, or (at your option) any later version. */

/* This library is distributed in the hope that it will be useful, */
/* but WITHOUT ANY WARRANTY; without even the implied warranty of */
/* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU */
/* Lesser General Public License for more details. */

/* You should have received a copy of the GNU Lesser General Public */
/* License along with this library; if not, write to the Free Software */
/* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA */

 // <h1>Example 14 - Laplace Equation in the L-Shaped Domain</h1>
 //
 // This example solves the Laplace equation on the classic "L-shaped"
 // domain with adaptive mesh refinement.  In this case, the exact
 // solution is u(r,\theta) = r^{2/3} * \sin ( (2/3) * \theta), but
 // the standard Kelly error indicator is used to estimate the error.
 // The initial mesh contains three QUAD9 elements which represent the
 // standard 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 file
 // named "lshaped.xda".  In addition, an input file named "ex14.in"
 // is provided which allows the user to set several parameters for
 // the solution so that the problem can be re-run without a
 // re-compile.  The solution technique employed is to have a
 // refinement loop with a linear solve inside followed by a
 // refinement of the grid and projection of the solution to the new grid
 // In the final loop iteration, there is no additional
 // refinement 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 of
 // elements which will be refined / coarsened at each step.

// LibMesh include files.
#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"

// Function prototype.  This is the function that will assemble
// the linear system for our Laplace problem.  Note that the
// function will take the \p EquationSystems object and the
// name of the system we are assembling as input.  From the
// \p EquationSystems object we have acess to the \p Mesh and
// other objects we might need.
void assemble_laplace(EquationSystems& es,
                      const std::string& system_name);


// Prototype for calculation of the exact solution.  Useful
// for setting boundary conditions.
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);

// Prototype for calculation of the gradient of the exact solution.  
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);


// These are non-const because the input file may change it,
// It is global because our exact_* functions use it.

// Set the dimensionality of the mesh
unsigned int dim = 2;

// Choose whether or not to use the singular solution
bool singularity = true;


int main(int argc, char** argv)
{
  // Initialize libMesh.
  LibMeshInit init (argc, argv);

#ifndef ENABLE_AMR
  if (libMesh::processor_id() == 0)
    std::cerr << "ERROR: This example requires libMesh to be\n"
              << "compiled with AMR support!"
              << std::endl;
  return 0;
#else

  // Parse the input file
  GetPot input_file("ex14.in");

  // Read in parameters from the input file
  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);
  
  // Output file for plotting the error as a function of
  // the number of degrees of freedom.
  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;
  
  // Create an n-dimensional mesh.
  Mesh mesh (dim);
  
  // Read in the mesh
  if (dim == 1)
    MeshTools::Generation::build_line(mesh,1,-1.,0.);
  else if (dim == 2)
    mesh.read("lshaped.xda");
  else
    mesh.read("lshaped3D.xda");

  // Use triangles if the config file says so
  if (element_type == "simplex")
    MeshTools::Modification::all_tri(mesh);

  // We used first order elements to describe the geometry,
  // but we may need second order elements to hold the degrees
  // of freedom
  if (approx_order > 1 || refine_type != "h")
    mesh.all_second_order();

  // Mesh Refinement object
  MeshRefinement mesh_refinement(mesh);
  mesh_refinement.refine_fraction() = refine_percentage;
  mesh_refinement.coarsen_fraction() = coarsen_percentage;
  mesh_refinement.max_h_level() = max_r_level;

  // Create an equation systems object.
  EquationSystems equation_systems (mesh);

  // Declare the system and its variables.
  // Creates a system named "Laplace"
  LinearImplicitSystem& system =
    equation_systems.add_system<LinearImplicitSystem> ("Laplace");
  
  // Adds the variable "u" to "Laplace", using 
  // the finite element type and order specified
  // in the config file
  system.add_variable("u", static_cast<Order>(approx_order),
                      Utility::string_to_enum<FEFamily>(approx_type));

  // Give the system a pointer to the matrix assembly
  // function.
  system.attach_assemble_function (assemble_laplace);

  // Initialize the data structures for the equation system.
  equation_systems.init();

  // Set linear solver max iterations
  equation_systems.parameters.set<unsigned int>("linear solver maximum iterations")
    = max_linear_iterations;

  // Linear solver tolerance.
  equation_systems.parameters.set<Real>("linear solver tolerance") =
    TOLERANCE * TOLERANCE * TOLERANCE;
  
  // Prints information about the system to the screen.
  equation_systems.print_info();

  // Construct ExactSolution object and attach solution functions
  ExactSolution exact_sol(equation_systems);
  exact_sol.attach_exact_value(exact_solution);
  exact_sol.attach_exact_deriv(exact_derivative);

  // Use higher quadrature order for more accurate error results
  exact_sol.extra_quadrature_order(extra_error_quadrature);

  // A refinement loop.
  for (unsigned int r_step=0; r_step<max_r_steps; r_step++)
    {
      std::cout << "Beginning Solve " << r_step << std::endl;
      
      // Solve the system "Laplace", just like example 2.
      system.solve();

      std::cout << "System has: " << equation_systems.n_active_dofs()
                << " degrees of freedom."
                << std::endl;

      std::cout << "Linear solver converged at step: "
                << system.n_linear_iterations()
                << ", final residual: "
                << system.final_linear_residual()
                << std::endl;
      
      // Compute the error.
      exact_sol.compute_error("Laplace", "u");

      // Print out the error values
      std::cout << "L2-Error is: "
                << exact_sol.l2_error("Laplace", "u")
                << std::endl;
      std::cout << "H1-Error is: "
                << exact_sol.h1_error("Laplace", "u")
                << std::endl;

      // Print to output file
      out << equation_systems.n_active_dofs() << " "
          << exact_sol.l2_error("Laplace", "u") << " "
          << exact_sol.h1_error("Laplace", "u") << std::endl;

      // Possibly refine the mesh
      if (r_step+1 != max_r_steps)
        {
          std::cout << "  Refining the mesh..." << std::endl;

          if (uniform_refine == 0)
            {

              // The \p ErrorVector is a particular \p StatisticsVector
              // for computing error information on a finite element mesh.
              ErrorVector error;
              
              if (indicator_type == "exact")
                {
                  // 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.
                  ExactErrorEstimator error_estimator;

                  error_estimator.attach_exact_value(exact_solution);
                  error_estimator.attach_exact_deriv(exact_derivative);

                  // We optimize in H1 norm
                  error_estimator.sobolev_order() = 1;

                  // 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.
                  error_estimator.estimate_error (system, error);
                }
              else if (indicator_type == "patch")
                {
                  // The patch recovery estimator should give a
                  // good estimate of the solution interpolation
                  // error.
                  PatchRecoveryErrorEstimator error_estimator;

                  error_estimator.estimate_error (system, error);
                }
              else if (indicator_type == "uniform")
                {
                  // Error indication based on uniform refinement
                  // is reliable, but very expensive.
                  UniformRefinementEstimator error_estimator;

                  error_estimator.estimate_error (system, error);
                }
              else
                {
                  libmesh_assert (indicator_type == "kelly");

                  // 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
                  KellyErrorEstimator error_estimator;

                  error_estimator.estimate_error (system, error);
                }
              
              // 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.
              mesh_refinement.flag_elements_by_error_fraction (error);

              // If we are doing adaptive p refinement, we want
              // elements flagged for that instead.
              if (refine_type == "p")
                mesh_refinement.switch_h_to_p_refinement();
              // If we are doing "matched hp" refinement, we
              // flag elements for both h and p
              if (refine_type == "matchedhp")
                mesh_refinement.add_p_to_h_refinement();
              // If we are doing hp refinement, we 
              // try switching some elements from h to p
              if (refine_type == "hp")
                {
                  HPCoarsenTest hpselector;
                  hpselector.select_refinement(system);
                }
              // If we are doing "singular hp" refinement, we 
              // try switching most elements from h to p
              if (refine_type == "singularhp")
                {
                  // This only differs from p refinement for
                  // the singular problem
                  libmesh_assert (singularity);
                  HPSingularity hpselector;
                  // Our only singular point is at the origin
                  hpselector.singular_points.push_back(Point());
                  hpselector.select_refinement(system);
                }
              
              // This call actually refines and coarsens the flagged
              // elements.
              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);