/* $Id: ex11.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 11 - Stokes Equations - Systems of Equations</h1> // // This example shows how a simple, linear system of equations // can be solved in parallel.  The system of equations are the familiar // Stokes equations for low-speed incompressible fluid flow.// C++ include files that we need#include <iostream>#include <algorithm>#include <math.h>// Basic include file needed for the mesh functionality.#include "libmesh.h"#include "mesh.h"#include "mesh_generation.h"#include "gmv_io.h"#include "equation_systems.h"#include "fe.h"#include "quadrature_gauss.h"#include "dof_map.h"#include "sparse_matrix.h"#include "numeric_vector.h"#include "dense_matrix.h"#include "dense_vector.h"#include "linear_implicit_system.h"// For systems of equations the \p DenseSubMatrix// and \p DenseSubVector provide convenient ways for// assembling the element matrix and vector on a// component-by-component basis.#include "dense_submatrix.h"#include "dense_subvector.h"// The definition of a geometric element#include "elem.h"// Function prototype.  This function will assemble the system// matrix and right-hand-side.void assemble_stokes (EquationSystems& es,                      const std::string& system_name);// The main program.int main (int argc, char** argv){  // Initialize libMesh.  LibMeshInit init (argc, argv);  // Set the dimensionality of the mesh = 2  const unsigned int dim = 2;           // Create a two-dimensional mesh.  Mesh mesh (dim);      // Use the MeshTools::Generation mesh generator to create a uniform  // grid on the square [-1,1]^D.  We instruct the mesh generator  // to build a mesh of 8x8 \p Quad9 elements in 2D, or \p Hex27  // elements in 3D.  Building these higher-order elements allows  // us to use higher-order approximation, as in example 3.  MeshTools::Generation::build_square (mesh,                                       15, 15,                                       0., 1.,                                       0., 1.,                                       QUAD9);    // Print information about the mesh to the screen.  mesh.print_info();    // Create an equation systems object.  EquationSystems equation_systems (mesh);    // Declare the system and its variables.  // Create a transient system named "Convection-Diffusion"  LinearImplicitSystem & system =     equation_systems.add_system<LinearImplicitSystem> ("Stokes");    // Add the variables "u" & "v" to "Stokes".  They  // will be approximated using second-order approximation.  system.add_variable ("u", SECOND);  system.add_variable ("v", SECOND);  // Add the variable "p" to "Stokes". This will  // be approximated with a first-order basis,  // providing an LBB-stable pressure-velocity pair.  system.add_variable ("p", FIRST);  // Give the system a pointer to the matrix assembly  // function.  system.attach_assemble_function (assemble_stokes);    // Initialize the data structures for the equation system.  equation_systems.init ();  equation_systems.parameters.set<unsigned int>("linear solver maximum iterations") = 250;  equation_systems.parameters.set<Real>        ("linear solver tolerance") = TOLERANCE;        // Prints information about the system to the screen.  equation_systems.print_info();      // Assemble & solve the linear system,  // then write the solution.  equation_systems.get_system("Stokes").solve();  GMVIO(mesh).write_equation_systems ("out.gmv",                                      equation_systems);  // All done.    return 0;}void assemble_stokes (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 == "Stokes");    // 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 Convection-Diffusion system object.  LinearImplicitSystem & system =    es.get_system<LinearImplicitSystem> ("Stokes");  // Numeric ids corresponding to each variable in the system  const unsigned int u_var = system.variable_number ("u");  const unsigned int v_var = system.variable_number ("v");  const unsigned int p_var = system.variable_number ("p");    // Get the Finite Element type for "u".  Note this will be  // the same as the type for "v".  FEType fe_vel_type = system.variable_type(u_var);    // Get the Finite Element type for "p".  FEType fe_pres_type = system.variable_type(p_var);  // Build a Finite Element object of the specified type for  // the velocity variables.  AutoPtr<FEBase> fe_vel  (FEBase::build(dim, fe_vel_type));      // Build a Finite Element object of the specified type for  // the pressure variables.  AutoPtr<FEBase> fe_pres (FEBase::build(dim, fe_pres_type));    // A Gauss quadrature rule for numerical integration.  // Let the \p FEType object decide what order rule is appropriate.  QGauss qrule (dim, fe_vel_type.default_quadrature_order());  // Tell the finite element objects to use our quadrature rule.  fe_vel->attach_quadrature_rule (&qrule);  fe_pres->attach_quadrature_rule (&qrule);    // 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_vel->get_JxW();    // The element shape function gradients for the velocity  // variables evaluated at the quadrature points.  const std::vector<std::vector<RealGradient> >& dphi = fe_vel->get_dphi();  // The element shape functions for the pressure variable  // evaluated at the quadrature points.  const std::vector<std::vector<Real> >& psi = fe_pres->get_phi();    // A reference to the \p DofMap object for this system.  The \p DofMap  // object handles the index translation from node and element numbers  // to degree of freedom numbers.  We will talk more about the \p DofMap  // in future examples.  const DofMap & dof_map = system.get_dof_map();  // Define data structures to contain the element matrix  // and right-hand-side vector contribution.  Following  // basic finite element terminology we will denote these  // "Ke" and "Fe".  DenseMatrix<Number> Ke;  DenseVector<Number> Fe;  DenseSubMatrix<Number>