/* $Id: ex3.C 2837 2008-05-08 17:23:37Z roystgnr $ *//* The Next Great Finite Element Library. *//* Copyright (C) 2003  Benjamin S. Kirk *//* 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 3 - Solving a Poisson Problem</h1> // // This is the third example program.  It builds on // the second example program by showing how to solve a simple // Poisson system.  This example also introduces the notion // of customized matrix assembly functions, working with an // exact solution, and using element iterators. // We will not comment on things that // were already explained in the second example.// C++ include files that we need#include <iostream>#include <algorithm>#include <math.h>// Basic include files needed for the mesh functionality.#include "libmesh.h"#include "mesh.h"#include "mesh_generation.h"#include "gmv_io.h"#include "linear_implicit_system.h"#include "equation_systems.h"// Define the Finite Element object.#include "fe.h"// Define Gauss quadrature rules.#include "quadrature_gauss.h"// Define useful datatypes for finite element// matrix and vector components.#include "sparse_matrix.h"#include "numeric_vector.h"#include "dense_matrix.h"#include "dense_vector.h"#include "elem.h"// Define the DofMap, which handles degree of freedom// indexing.#include "dof_map.h"// Function prototype.  This is the function that will assemble// the linear system for our Poisson problem.  Note that the// function will take the  EquationSystems object and the// name of the system we are assembling as input.  From the//  EquationSystems object we have access to the  Mesh and// other objects we might need.void assemble_poisson(EquationSystems& es,                      const std::string& system_name);// Function prototype for the exact solution.Real exact_solution (const Real x,                     const Real y,                     const Real z = 0.);int main (int argc, char** argv){  // Initialize libraries, like in example 2.  LibMeshInit init (argc, argv);  // Brief message to the user regarding the program name  // and command line arguments.  std::cout << "Running " << argv[0];    for (int i=1; i<argc; i++)    std::cout << " " << argv[i];    std::cout << std::endl << std::endl;    // Create a 2D mesh.  Mesh mesh (2);      // Use the MeshTools::Generation mesh generator to create a uniform  // grid on the square [-1,1]^2.  We instruct the mesh generator  // to build a mesh of 15x15 QUAD9 elements.  Building QUAD9  // elements instead of the default QUAD4's we used in example 2  // allow us to use higher-order approximation.  MeshTools::Generation::build_square (mesh,                                        15, 15,                                       -1., 1.,                                       -1., 1.,                                       QUAD9);  // Print information about the mesh to the screen.  // Note that 5x5 QUAD9 elements actually has 11x11 nodes,  // so this mesh is significantly larger than the one in example 2.  mesh.print_info();    // Create an equation systems object.  EquationSystems equation_systems (mesh);    // Declare the Poisson system and its variables.  // The Poisson system is another example of a steady system.  equation_systems.add_system<LinearImplicitSystem> ("Poisson");  // Adds the variable "u" to "Poisson".  "u"  // will be approximated using second-order approximation.  equation_systems.get_system("Poisson").add_variable("u", SECOND);  // Give the system a pointer to the matrix assembly  // function.  This will be called when needed by the  // library.  equation_systems.get_system("Poisson").attach_assemble_function (assemble_poisson);    // Initialize the data structures for the equation system.  equation_systems.init();    // Prints information about the system to the screen.  equation_systems.print_info();  // Solve the system "Poisson".  Note that calling this  // member will assemble the linear system and invoke  // the default Petsc solver, however the solver can be  // controlled from the command line.  For example,  // you can invoke conjugate gradient with:  //  // ./ex3 -ksp_type cg  //  // and you can get a nice X-window that monitors the solver  // convergence with:  //  // ./ex3 -ksp_xmonitor  //  // if you linked against the appropriate X libraries when you  // built PETSc.  equation_systems.get_system("Poisson").solve();  // After solving the system write the solution  // to a GMV-formatted plot file.  GMVIO (mesh).write_equation_systems ("out.gmv", equation_systems);  // All done.    return 0;}// We now define the matrix assembly function for the// Poisson system.  We need to first compute element// matrices and right-hand sides, and then take into// account the boundary conditions, which will be handled// via a penalty method.void assemble_poisson(EquationSystems& es,                      const std::string& system_name){    // It is a good idea to make sure we are assembling  // the proper system.  libmesh_assert (system_name == "Poisson");    // Get a constant reference to the mesh object.  const MeshBase& mesh = es.get_mesh();  // The dimension that we are running  const unsigned int dim = mesh.mesh_dimension();  // Get a reference to the LinearImplicitSystem we are solving  LinearImplicitSystem& system = es.get_system<LinearImplicitSystem> ("Poisson");  // A reference to the  DofMap object for this system.  The  DofMap  // object handles the index translation from node and element numbers  // to degree of freedom numbers.  We will talk more about the  DofMap  // in future examples.  const DofMap& dof_map = system.get_dof_map();    // Get a constant reference to the Finite Element type  // for the first (and only) variable in the system.  FEType fe_type = dof_map.variable_type(0);    // Build a Finite Element object of the specified type.  Since the  // FEBase::build() member dynamically creates memory we will  // store the object as an AutoPtr<FEBase>.  This can be thought  // of as a pointer that will clean up after itself.  Example 4  // describes some advantages of  AutoPtr's in the context of  // quadrature rules.  AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));    // A 5th order Gauss quadrature rule for numerical integration.  QGauss qrule (dim, FIFTH);    // Tell the finite element object to use our quadrature rule.  fe->attach_quadrature_rule (&qrule);    // Declare a special finite element object for  // boundary integration.  AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));    // Boundary integration requires one quadraure rule,  // with dimensionality one less than the dimensionality  // of the element.  QGauss qface(dim-1, FIFTH);    // Tell the finite element object to use our  // quadrature rule.  fe_face->attach_quadrature_rule (&qface);  // Here we define some references to cell-specific data that  // will be used to assemble the linear system.  //  // The element Jacobian * quadrature weight at each integration point.     const std::vector<Real>& JxW = fe->get_JxW();  // The physical XY locations of the quadrature points on the element.  // These might be useful for evaluating spatially varying material