ex15.c

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

/* $Id: ex15.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 15 - Biharmonic Equation</h1>
 //
 // This example solves the Biharmonic equation on a square or cube,
 // using a Galerkin formulation with C1 elements approximating the
 // H^2_0 function space.
 // The initial mesh contains two TRI6, one QUAD9 or one HEX27
 // An input file named "ex15.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 "ex15.in", the variable
 // "max_r_steps" controls the number of refinement steps, and
 // "max_r_level" controls the maximum element refinement level.

// LibMesh include files.
#include "mesh.h"
#include "equation_systems.h"
#include "linear_implicit_system.h"
#include "gmv_io.h"
#include "fe.h"
#include "quadrature.h"
#include "dense_matrix.h"
#include "dense_vector.h"
#include "sparse_matrix.h"
#include "mesh_generation.h"
#include "mesh_modification.h"
#include "mesh_refinement.h"
#include "error_vector.h"
#include "fourth_error_estimators.h"
#include "getpot.h"
#include "exact_solution.h"
#include "dof_map.h"
#include "numeric_vector.h"
#include "elem.h"
#include "tensor_value.h"

// Function prototype.  This is the function that will assemble
// the linear system for our Biharmonic 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_biharmonic(EquationSystems& es,
                      const std::string& system_name);


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

Number exact_3D_solution(const Point& p,
  const Parameters&, const std::string&, const std::string&);

// Prototypes for calculation of the gradient of the exact solution.  
// Necessary for setting boundary conditions in H^2_0 and testing
// H^1 convergence of the solution
Gradient exact_2D_derivative(const Point& p,
  const Parameters&, const std::string&, const std::string&);

Gradient exact_3D_derivative(const Point& p,
  const Parameters&, const std::string&, const std::string&);

Tensor exact_2D_hessian(const Point& p,
  const Parameters&, const std::string&, const std::string&);

Tensor exact_3D_hessian(const Point& p,
  const Parameters&, const std::string&, const std::string&);

Number forcing_function_2D(const Point& p);

Number forcing_function_3D(const Point& p);

// Pointers to dimension-independent functions
Number (*exact_solution)(const Point& p,
  const Parameters&, const std::string&, const std::string&);
Gradient (*exact_derivative)(const Point& p,
  const Parameters&, const std::string&, const std::string&);
Tensor (*exact_hessian)(const Point& p,
  const Parameters&, const std::string&, const std::string&);
Number (*forcing_function)(const Point& p);



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

#ifndef ENABLE_SECOND_DERIVATIVES
  if (libMesh::processor_id() == 0)
    std::cerr << "ERROR: This example requires the library to be "
              << "compiled with second derivatives support!"
              << std::endl;
  return 0;
#else

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

  // Read in parameters from the input file
  const unsigned int max_r_level = input_file("max_r_level", 10);
  const unsigned int max_r_steps = input_file("max_r_steps", 4);
  const std::string approx_type  = input_file("approx_type",
                                              "HERMITE");
  const unsigned int uniform_refine =
                  input_file("uniform_refine", 0);
  const Real refine_percentage =
                  input_file("refine_percentage", 0.5);
  const Real coarsen_percentage =
                  input_file("coarsen_percentage", 0.5);
  const unsigned int dim =
                  input_file("dimension", 2);
  const unsigned int max_linear_iterations =
                  input_file("max_linear_iterations", 10000);

  // We have only defined 2 and 3 dimensional problems
  libmesh_assert (dim == 2 || dim == 3);

  // Currently only the Hermite cubics give a 3D C^1 basis
  libmesh_assert (dim == 2 || approx_type == "HERMITE");

  // Create a dim-dimensional mesh.
  Mesh mesh (dim);
  
  // Output file for plotting the error 
  std::string output_file = "";

  if (dim == 2)
    output_file += "2D_";
  else if (dim == 3)
    output_file += "3D_";

  if (approx_type == "HERMITE")
    output_file += "hermite_";
  else if (approx_type == "SECOND")
    output_file += "reducedclough_";
  else
    output_file += "clough_";

  if (uniform_refine == 0)
    output_file += "adaptive";
  else
    output_file += "uniform";

  std::string gmv_file = output_file;
  gmv_file += ".gmv";
  output_file += ".m";

  std::ofstream out (output_file.c_str());
  out << "% dofs     L2-error     H1-error      H2-error\n"
      << "e = [\n";
  
  // Set up the dimension-dependent coarse mesh and solution
  // We build more than one cell so as to avoid bugs on fewer than 
  // 4 processors in 2D or 8 in 3D.
  if (dim == 2)
    {
      MeshTools::Generation::build_square(mesh, 2, 2);
      exact_solution = &exact_2D_solution;
      exact_derivative = &exact_2D_derivative;
      exact_hessian = &exact_2D_hessian;
      forcing_function = &forcing_function_2D;
    }
  else if (dim == 3)
    {
      MeshTools::Generation::build_cube(mesh, 2, 2, 2);
      exact_solution = &exact_3D_solution;
      exact_derivative = &exact_3D_derivative;
      exact_hessian = &exact_3D_hessian;
      forcing_function = &forcing_function_3D;
    }

  // Convert the mesh to second order: necessary for computing with
  // Clough-Tocher elements, useful for getting slightly less 
  // broken gmv output with Hermite elements
  mesh.all_second_order();

  // Convert it to triangles if necessary
  if (approx_type != "HERMITE")
    MeshTools::Modification::all_tri(mesh);

  // 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.
  // Create a system named "Biharmonic"
  LinearImplicitSystem& system =
    equation_systems.add_system<LinearImplicitSystem> ("Biharmonic");

  // Adds the variable "u" to "Biharmonic".  "u"
  // will be approximated using Hermite tensor product squares
  // or (possibly reduced) cubic Clough-Tocher triangles
  if (approx_type == "HERMITE")
    system.add_variable("u", THIRD, HERMITE);
  else if (approx_type == "SECOND")
    system.add_variable("u", SECOND, CLOUGH);
  else if (approx_type == "CLOUGH")
    system.add_variable("u", THIRD, CLOUGH);
  else
    libmesh_error();

  // Give the system a pointer to the matrix assembly
  // function.
  system.attach_assemble_function (assemble_biharmonic);
      
  // 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;
      
  // Prints information about the system to the screen.
  equation_systems.print_info();

  // Construct ExactSolution object and attach function to compute exact solution
  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);

  // Construct zero solution object, useful for computing solution norms
  // Attaching "zero_solution" functions is unnecessary
  ExactSolution zero_sol(equation_systems);

  // A refinement loop.
  for (unsigned int r_step=0; r_step<max_r_steps; r_step++)
    {
      mesh.print_info();
      equation_systems.print_info();

      std::cout << "Beginning Solve " << r_step << std::endl;
      
      // Solve the system "Biharmonic", just like example 2.
      system.solve();

      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("Biharmonic", "u");
      // Compute the norm.
      zero_sol.compute_error("Biharmonic", "u");

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

      // Print to output file
      out << equation_systems.n_active_dofs() << " "
          << exact_sol.l2_error("Biharmonic", "u") << " "
          << exact_sol.h1_error("Biharmonic", "u") << " "
          << exact_sol.h2_error("Biharmonic", "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)
            {
              ErrorVector error;
              LaplacianErrorEstimator error_estimator;

              error_estimator.estimate_error(system, error);
              mesh_refinement.flag_elements_by_elem_fraction (error);

              std::cout << "Mean Error: " << error.mean() <<
                              std::endl;
              std::cout << "Error Variance: " << error.variance() <<
                              std::endl;

              mesh_refinement.refine_and_coarsen_elements();
            }
          else
            {
              mesh_refinement.uniformly_refine(1);
            }
              
          // 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.
          equation_systems.reinit ();
        }
    }            
  
  // Write out the solution
  // After solving the system write the solution
  // to a GMV-formatted plot file.
  GMVIO (mesh).write_equation_systems (gmv_file,
                                           equation_systems);

  // Close up the output file.
  out << "];\n"
      << "hold on\n"
      << "loglog(e(:,1), e(:,2), 'bo-');\n"
      << "loglog(e(:,1), e(:,3), 'ro-');\n"
      << "loglog(e(:,1), e(:,4), 'go-');\n"
      << "xlabel('log(dofs)');\n"
      << "ylabel('log(error)');\n"
      << "title('C1 " << approx_type << " elements');\n"
      << "legend('L2-error', 'H1-error', 'H2-error');\n";
  
  // All done.  
  return 0;
#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
{
  // x and y coordinates in space
  const Real x = p(0);
  const Real y = p(1);

  // analytic solution value
  return 256.*(x-x*x)*(x-x*x)*(y-y*y)*(y-y*y);
}


// We now define the gradient of the exact solution
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
{
  // x and y coordinates in space
  const Real x = p(0);
  const Real y = p(1);

  // First derivatives to be returned.
  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;
}


// We now define the hessian of the exact solution
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