The source file miscellaneous_ex3.C with comments:

Miscellaneous Example 3 - 2D Laplace-Young Problem Using Nonlinear Solvers



This example shows how to use the NonlinearImplicitSystem class to efficiently solve nonlinear problems in parallel.

In nonlinear systems, we aim at finding x that satisfy R(x) = 0. In nonlinear finite element analysis, the residual is typically of the form R(x) = K(x)*x - f, with K(x) the system matrix and f the "right-hand-side". The NonlinearImplicitSystem class expects two callback functions to compute the residual R and its Jacobian for the Newton iterations. Here, we just approximate the true Jacobian by K(x).

You can turn on preconditining of the matrix free system using the jacobian by passing "-pre" on the command line. Currently this only work with Petsc so this isn't used by using "make run"

This example also runs with the experimental Trilinos NOX solvers by specifying the --use-trilinos command line argument.



C++ include files that we need
        #include <iostream>
        #include <algorithm>
        #include <cmath>
        
Various include files needed for the mesh & solver functionality.
        #include "libmesh/libmesh.h"
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_refinement.h"
        #include "libmesh/exodusII_io.h"
        #include "libmesh/equation_systems.h"
        #include "libmesh/fe.h"
        #include "libmesh/quadrature_gauss.h"
        #include "libmesh/dof_map.h"
        #include "libmesh/sparse_matrix.h"
        #include "libmesh/numeric_vector.h"
        #include "libmesh/dense_matrix.h"
        #include "libmesh/dense_vector.h"
        #include "libmesh/elem.h"
        #include "libmesh/string_to_enum.h"
        #include "libmesh/getpot.h"
        
The nonlinear solver and system we will be using
        #include "libmesh/nonlinear_solver.h"
        #include "libmesh/nonlinear_implicit_system.h"
        
Necessary for programmatically setting petsc options
        #ifdef LIBMESH_HAVE_PETSC
        #include <petsc.h>
        #endif
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
A reference to our equation system
        EquationSystems *_equation_system = NULL;
        
Let-s define the physical parameters of the equation
        const Real kappa = 1.;
        const Real sigma = 0.2;
        
        
This function computes the Jacobian K(x)
        void compute_jacobian (const NumericVector<Number>& soln,
                               SparseMatrix<Number>&  jacobian,
                               NonlinearImplicitSystem& /*sys*/)
        {
Get a reference to the equation system.
          EquationSystems &es = *_equation_system;
        
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 NonlinearImplicitSystem we are solving
          NonlinearImplicitSystem& system =
            es.get_system<NonlinearImplicitSystem>("Laplace-Young");
        
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();
        
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 \p FEBase::build() member dynamically creates memory we will store the object as an \p AutoPtr. This can be thought of as a pointer that will clean up after itself.
          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);
        
Here we define some references to cell-specific data that will be used to assemble the linear system. We begin with the element Jacobian * quadrature weight at each integration point.
          const std::vector<Real>& JxW = fe->get_JxW();
        
The element shape functions evaluated at the quadrature points.
          const std::vector<std::vector<Real> >& phi = fe->get_phi();
        
The element shape function gradients evaluated at the quadrature points.
          const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
        
Define data structures to contain the Jacobian element matrix. Following basic finite element terminology we will denote these "Ke". More detail is in example 3.
          DenseMatrix<Number> Ke;
        
This vector will hold the degree of freedom indices for the element. These define where in the global system the element degrees of freedom get mapped.
          std::vector<dof_id_type> dof_indices;
        
Now we will loop over all the active elements in the mesh which are local to this processor. We will compute the element Jacobian contribution.
          MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
          const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
        
          for ( ; el != end_el; ++el)
            {
Store a pointer to the element we are currently working on. This allows for nicer syntax later.
              const Elem* elem = *el;
        
Get the degree of freedom indices for the current element. These define where in the global matrix and right-hand-side this element will contribute to.
              dof_map.dof_indices (elem, dof_indices);
        
Compute the element-specific data for the current element. This involves computing the location of the quadrature points (q_point) and the shape functions (phi, dphi) for the current element.
              fe->reinit (elem);
        
Zero the element Jacobian before summing them. We use the resize member here because the number of degrees of freedom might have changed from the last element. Note that this will be the case if the element type is different (i.e. the last element was a triangle, now we are on a quadrilateral).
              Ke.resize (dof_indices.size(),
                         dof_indices.size());
        
Now we will build the element Jacobian. This involves a double loop to integrate the test funcions (i) against the trial functions (j). Note that the Jacobian depends on the current solution x, which we access using the soln vector.

              for (unsigned int qp=0; qp<qrule.n_points(); qp++)
                {
                  Gradient grad_u;
        
                  for (unsigned int i=0; i<phi.size(); i++)
                    grad_u += dphi[i][qp]*soln(dof_indices[i]);
        
                  const Number K = 1./std::sqrt(1. + grad_u*grad_u);
        
                  for (unsigned int i=0; i<phi.size(); i++)
                    for (unsigned int j=0; j<phi.size(); j++)
                      Ke(i,j) += JxW[qp]*(
                                          K*(dphi[i][qp]*dphi[j][qp]) +
                                          kappa*phi[i][qp]*phi[j][qp]
                                          );
                }
        
              dof_map.constrain_element_matrix (Ke, dof_indices);
        
Add the element matrix to the system Jacobian.
              jacobian.add_matrix (Ke, dof_indices);
            }
        
That's it.
        }
        
        
Here we compute the residual R(x) = K(x)*x - f. The current solution x is passed in the soln vector
        void compute_residual (const NumericVector<Number>& soln,
                               NumericVector<Number>& residual,
                               NonlinearImplicitSystem& /*sys*/)
        {
          EquationSystems &es = *_equation_system;
        
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();
          libmesh_assert_equal_to (dim, 2);
        
Get a reference to the NonlinearImplicitSystem we are solving
          NonlinearImplicitSystem& system =
            es.get_system<NonlinearImplicitSystem>("Laplace-Young");
        
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();
        
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 \p FEBase::build() member dynamically creates memory we will store the object as an \p AutoPtr. This can be thought of as a pointer that will clean up after itself.
          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 finte 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. We begin with the element Jacobian * quadrature weight at each integration point.
          const std::vector<Real>& JxW = fe->get_JxW();
        
The element shape functions evaluated at the quadrature points.
          const std::vector<std::vector<Real> >& phi = fe->get_phi();
        
The element shape function gradients evaluated at the quadrature points.
          const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
        
Define data structures to contain the resdual contributions
          DenseVector<Number> Re;
        
This vector will hold the degree of freedom indices for the element. These define where in the global system the element degrees of freedom get mapped.
          std::vector<dof_id_type> dof_indices;
        
Now we will loop over all the active elements in the mesh which are local to this processor. We will compute the element residual.
          residual.zero();
        
          MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
          const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
        
          for ( ; el != end_el; ++el)
            {
Store a pointer to the element we are currently working on. This allows for nicer syntax later.
              const Elem* elem = *el;
        
Get the degree of freedom indices for the current element. These define where in the global matrix and right-hand-side this element will contribute to.
              dof_map.dof_indices (elem, dof_indices);
        
Compute the element-specific data for the current element. This involves computing the location of the quadrature points (q_point) and the shape functions (phi, dphi) for the current element.
              fe->reinit (elem);
        
We use the resize member here because the number of degrees of freedom might have changed from the last element. Note that this will be the case if the element type is different (i.e. the last element was a triangle, now we are on a quadrilateral).
              Re.resize (dof_indices.size());
        
Now we will build the residual. This involves the construction of the matrix K and multiplication of it with the current solution x. We rearrange this into two loops: In the first, we calculate only the contribution of K_ij*x_j which is independent of the row i. In the second loops, we multiply with the row-dependent part and add it to the element residual.

              for (unsigned int qp=0; qp<qrule.n_points(); qp++)
                {
                  Number u = 0;
                  Gradient grad_u;
        
                  for (unsigned int j=0; j<phi.size(); j++)
                    {
                      u      += phi[j][qp]*soln(dof_indices[j]);
                      grad_u += dphi[j][qp]*soln(dof_indices[j]);
                    }
        
                  const Number K = 1./std::sqrt(1. + grad_u*grad_u);
        
                  for (unsigned int i=0; i<phi.size(); i++)
                    Re(i) += JxW[qp]*(
                                      K*(dphi[i][qp]*grad_u) +
                                      kappa*phi[i][qp]*u
                                      );
                }
        
At this point the interior element integration has been completed. However, we have not yet addressed boundary conditions.

The following loops over the sides of the element. If the element has no neighbor on a side then that side MUST live on a boundary of the domain.
              for (unsigned int side=0; side<elem->n_sides(); side++)
                if (elem->neighbor(side) == NULL)
                  {
The value of the shape functions at the quadrature points.
                    const std::vector<std::vector<Real> >&  phi_face = fe_face->get_phi();
        
The Jacobian * Quadrature Weight at the quadrature points on the face.
                    const std::vector<Real>& JxW_face = fe_face->get_JxW();
        
Compute the shape function values on the element face.
                    fe_face->reinit(elem, side);
        
Loop over the face quadrature points for integration.
                    for (unsigned int qp=0; qp<qface.n_points(); qp++)
                      {
This is the right-hand-side contribution (f), which has to be subtracted from the current residual
                        for (unsigned int i=0; i<phi_face.size(); i++)
                          Re(i) -= JxW_face[qp]*sigma*phi_face[i][qp];
                      }
                  }
        
              dof_map.constrain_element_vector (Re, dof_indices);
              residual.add_vector (Re, dof_indices);
            }
        
That's it.
        }
        
        
        
Begin the main program.
        int main (int argc, char** argv)
        {
Initialize libMesh and any dependent libaries, like in example 2.
          LibMeshInit init (argc, argv);
        
        #if !defined(LIBMESH_HAVE_PETSC) && !defined(LIBMESH_HAVE_TRILINOS)
          if (libMesh::processor_id() == 0)
            std::cerr << "ERROR: This example requires libMesh to be\n"
                      << "compiled with nonlinear solver support from\n"
                      << "PETSc or Trilinos!"
                      << std::endl;
          return 0;
        #endif
        
        #ifndef LIBMESH_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
        
Create a GetPot object to parse the command line
          GetPot command_line (argc, argv);
        
Check for proper calling arguments.
          if (argc < 3)
            {
              if (libMesh::processor_id() == 0)
                std::cerr << "Usage:\n"
                          <<"\t " << argv[0] << " -r 2"
                          << std::endl;
        
This handy function will print the file name, line number, and then abort.
              libmesh_error();
            }
        
Brief message to the user regarding the program name and command line arguments.
          else
            {
              std::cout << "Running " << argv[0];
        
              for (int i=1; i<argc; i++)
                std::cout << " " << argv[i];
        
              std::cout << std::endl << std::endl;
            }
        
        
Read number of refinements
          int nr = 2;
          if ( command_line.search(1, "-r") )
            nr = command_line.next(nr);
        
Read FE order from command line
          std::string order = "FIRST";
          if ( command_line.search(2, "-Order", "-o") )
            order = command_line.next(order);
        
Read FE Family from command line
          std::string family = "LAGRANGE";
          if ( command_line.search(2, "-FEFamily", "-f") )
            family = command_line.next(family);
        
Cannot use dicontinuous basis.
          if ((family == "MONOMIAL") || (family == "XYZ"))
            {
              std::cout << "ex19 currently requires a C^0 (or higher) FE basis." << std::endl;
              libmesh_error();
            }
        
          if ( command_line.search(1, "-pre") )
            {
        #ifdef LIBMESH_HAVE_PETSC
Use the jacobian for preconditioning.
              PetscOptionsSetValue("-snes_mf_operator",PETSC_NULL);
        #else
              std::cerr<<"Must be using PetsC to use jacobian based preconditioning"<<std::endl;
        
returning zero so that "make run" won't fail if we ever enable this capability there.
              return 0;
        #endif //LIBMESH_HAVE_PETSC
            }
        
Skip this 2D example if libMesh was compiled as 1D-only.
          libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
        
Create a mesh, with dimension to be overridden by the file, distributed across the default MPI communicator.
          Mesh mesh(init.comm());
        
          mesh.read ("lshaped.xda");
        
          if (order != "FIRST")
            mesh.all_second_order();
        
          MeshRefinement(mesh).uniformly_refine(nr);
        
Print information about the mesh to the screen.
          mesh.print_info();
        
Create an equation systems object.
          EquationSystems equation_systems (mesh);
          _equation_system = &equation_systems;
        
Declare the system and its variables.

Creates a system named "Laplace-Young"
          NonlinearImplicitSystem& system =
            equation_systems.add_system<NonlinearImplicitSystem> ("Laplace-Young");
        
        
Here we specify the tolerance for the nonlinear solver and the maximum of nonlinear iterations.
          equation_systems.parameters.set<Real>         ("nonlinear solver tolerance")          = 1.e-12;
          equation_systems.parameters.set<unsigned int> ("nonlinear solver maximum iterations") = 50;
        
        
Adds the variable "u" to "Laplace-Young". "u" will be approximated using second-order approximation.
          system.add_variable("u",
        		      Utility::string_to_enum<Order>   (order),
        		      Utility::string_to_enum<FEFamily>(family));
        
Give the system a pointer to the functions that update the residual and Jacobian.
          system.nonlinear_solver->residual = compute_residual;
          system.nonlinear_solver->jacobian = compute_jacobian;
        
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 "Laplace-Young", print the number of iterations and final residual
          equation_systems.get_system("Laplace-Young").solve();
        
Print out final convergence information. This duplicates some output from during the solve itself, but demonstrates another way to get this information after the solve is complete.
          std::cout << "Laplace-Young system solved at nonlinear iteration "
        	    << system.n_nonlinear_iterations()
        	    << " , final nonlinear residual norm: "
        	    << system.final_nonlinear_residual()
        	    << std::endl;
        
        #ifdef LIBMESH_HAVE_EXODUS_API
After solving the system write the solution
          ExodusII_IO (mesh).write_equation_systems ("out.e",
                                               equation_systems);
        #endif // #ifdef LIBMESH_HAVE_EXODUS_API
        #endif // #ifndef LIBMESH_ENABLE_AMR
        
All done.
          return 0;
        }



The source file miscellaneous_ex3.C without comments:

 
  
  
  #include <iostream>
  #include <algorithm>
  #include <cmath>
  
  #include "libmesh/libmesh.h"
  #include "libmesh/mesh.h"
  #include "libmesh/mesh_refinement.h"
  #include "libmesh/exodusII_io.h"
  #include "libmesh/equation_systems.h"
  #include "libmesh/fe.h"
  #include "libmesh/quadrature_gauss.h"
  #include "libmesh/dof_map.h"
  #include "libmesh/sparse_matrix.h"
  #include "libmesh/numeric_vector.h"
  #include "libmesh/dense_matrix.h"
  #include "libmesh/dense_vector.h"
  #include "libmesh/elem.h"
  #include "libmesh/string_to_enum.h"
  #include "libmesh/getpot.h"
  
  #include "libmesh/nonlinear_solver.h"
  #include "libmesh/nonlinear_implicit_system.h"
  
  #ifdef LIBMESH_HAVE_PETSC
  #include <petsc.h>
  #endif
  
  using namespace libMesh;
  
  EquationSystems *_equation_system = NULL;
  
  const Real kappa = 1.;
  const Real sigma = 0.2;
  
  
  void compute_jacobian (const NumericVector<Number>& soln,
                         SparseMatrix<Number>&  jacobian,
                         NonlinearImplicitSystem& /*sys*/)
  {
    EquationSystems &es = *_equation_system;
  
    const MeshBase& mesh = es.get_mesh();
  
    const unsigned int dim = mesh.mesh_dimension();
  
    NonlinearImplicitSystem& system =
      es.get_system<NonlinearImplicitSystem>("Laplace-Young");
  
    const DofMap& dof_map = system.get_dof_map();
  
    FEType fe_type = dof_map.variable_type(0);
  
    AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
  
    QGauss qrule (dim, FIFTH);
  
    fe->attach_quadrature_rule (&qrule);
  
    const std::vector<Real>& JxW = fe->get_JxW();
  
    const std::vector<std::vector<Real> >& phi = fe->get_phi();
  
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
    DenseMatrix<Number> Ke;
  
    std::vector<dof_id_type> dof_indices;
  
    MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
  
    for ( ; el != end_el; ++el)
      {
        const Elem* elem = *el;
  
        dof_map.dof_indices (elem, dof_indices);
  
        fe->reinit (elem);
  
        Ke.resize (dof_indices.size(),
                   dof_indices.size());
  
        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
          {
            Gradient grad_u;
  
            for (unsigned int i=0; i<phi.size(); i++)
              grad_u += dphi[i][qp]*soln(dof_indices[i]);
  
            const Number K = 1./std::sqrt(1. + grad_u*grad_u);
  
            for (unsigned int i=0; i<phi.size(); i++)
              for (unsigned int j=0; j<phi.size(); j++)
                Ke(i,j) += JxW[qp]*(
                                    K*(dphi[i][qp]*dphi[j][qp]) +
                                    kappa*phi[i][qp]*phi[j][qp]
                                    );
          }
  
        dof_map.constrain_element_matrix (Ke, dof_indices);
  
        jacobian.add_matrix (Ke, dof_indices);
      }
  
  }
  
  
  void compute_residual (const NumericVector<Number>& soln,
                         NumericVector<Number>& residual,
                         NonlinearImplicitSystem& /*sys*/)
  {
    EquationSystems &es = *_equation_system;
  
    const MeshBase& mesh = es.get_mesh();
  
    const unsigned int dim = mesh.mesh_dimension();
    libmesh_assert_equal_to (dim, 2);
  
    NonlinearImplicitSystem& system =
      es.get_system<NonlinearImplicitSystem>("Laplace-Young");
  
    const DofMap& dof_map = system.get_dof_map();
  
    FEType fe_type = dof_map.variable_type(0);
  
    AutoPtr<FEBase> fe (FEBase::build(dim, fe_type));
  
    QGauss qrule (dim, FIFTH);
  
    fe->attach_quadrature_rule (&qrule);
  
    AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type));
  
    QGauss qface(dim-1, FIFTH);
  
    fe_face->attach_quadrature_rule (&qface);
  
    const std::vector<Real>& JxW = fe->get_JxW();
  
    const std::vector<std::vector<Real> >& phi = fe->get_phi();
  
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
    DenseVector<Number> Re;
  
    std::vector<dof_id_type> dof_indices;
  
    residual.zero();
  
    MeshBase::const_element_iterator       el     = mesh.active_local_elements_begin();
    const MeshBase::const_element_iterator end_el = mesh.active_local_elements_end();
  
    for ( ; el != end_el; ++el)
      {
        const Elem* elem = *el;
  
        dof_map.dof_indices (elem, dof_indices);
  
        fe->reinit (elem);
  
        Re.resize (dof_indices.size());
  
  
        for (unsigned int qp=0; qp<qrule.n_points(); qp++)
          {
            Number u = 0;
            Gradient grad_u;
  
            for (unsigned int j=0; j<phi.size(); j++)
              {
                u      += phi[j][qp]*soln(dof_indices[j]);
                grad_u += dphi[j][qp]*soln(dof_indices[j]);
              }
  
            const Number K = 1./std::sqrt(1. + grad_u*grad_u);
  
            for (unsigned int i=0; i<phi.size(); i++)
              Re(i) += JxW[qp]*(
                                K*(dphi[i][qp]*grad_u) +
                                kappa*phi[i][qp]*u
                                );
          }
  
  
        for (unsigned int side=0; side<elem->n_sides(); side++)
          if (elem->neighbor(side) == NULL)
            {
              const std::vector<std::vector<Real> >&  phi_face = fe_face->get_phi();
  
              const std::vector<Real>& JxW_face = fe_face->get_JxW();
  
              fe_face->reinit(elem, side);
  
              for (unsigned int qp=0; qp<qface.n_points(); qp++)
                {
                  for (unsigned int i=0; i<phi_face.size(); i++)
                    Re(i) -= JxW_face[qp]*sigma*phi_face[i][qp];
                }
            }
  
        dof_map.constrain_element_vector (Re, dof_indices);
        residual.add_vector (Re, dof_indices);
      }
  
  }
  
  
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
  
  #if !defined(LIBMESH_HAVE_PETSC) && !defined(LIBMESH_HAVE_TRILINOS)
    if (libMesh::processor_id() == 0)
      std::cerr << "ERROR: This example requires libMesh to be\n"
                << "compiled with nonlinear solver support from\n"
                << "PETSc or Trilinos!"
                << std::endl;
    return 0;
  #endif
  
  #ifndef LIBMESH_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
  
    GetPot command_line (argc, argv);
  
    if (argc < 3)
      {
        if (libMesh::processor_id() == 0)
          std::cerr << "Usage:\n"
                    <<"\t " << argv[0] << " -r 2"
                    << std::endl;
  
        libmesh_error();
      }
  
    else
      {
        std::cout << "Running " << argv[0];
  
        for (int i=1; i<argc; i++)
          std::cout << " " << argv[i];
  
        std::cout << std::endl << std::endl;
      }
  
  
    int nr = 2;
    if ( command_line.search(1, "-r") )
      nr = command_line.next(nr);
  
    std::string order = "FIRST";
    if ( command_line.search(2, "-Order", "-o") )
      order = command_line.next(order);
  
    std::string family = "LAGRANGE";
    if ( command_line.search(2, "-FEFamily", "-f") )
      family = command_line.next(family);
  
    if ((family == "MONOMIAL") || (family == "XYZ"))
      {
        std::cout << "ex19 currently requires a C^0 (or higher) FE basis." << std::endl;
        libmesh_error();
      }
  
    if ( command_line.search(1, "-pre") )
      {
  #ifdef LIBMESH_HAVE_PETSC
        PetscOptionsSetValue("-snes_mf_operator",PETSC_NULL);
  #else
        std::cerr<<"Must be using PetsC to use jacobian based preconditioning"<<std::endl;
  
        return 0;
  #endif //LIBMESH_HAVE_PETSC
      }
  
    libmesh_example_assert(2 <= LIBMESH_DIM, "2D support");
  
    Mesh mesh(init.comm());
  
    mesh.read ("lshaped.xda");
  
    if (order != "FIRST")
      mesh.all_second_order();
  
    MeshRefinement(mesh).uniformly_refine(nr);
  
    mesh.print_info();
  
    EquationSystems equation_systems (mesh);
    _equation_system = &equation_systems;
  
  
    NonlinearImplicitSystem& system =
      equation_systems.add_system<NonlinearImplicitSystem> ("Laplace-Young");
  
  
    equation_systems.parameters.set<Real>         ("nonlinear solver tolerance")          = 1.e-12;
    equation_systems.parameters.set<unsigned int> ("nonlinear solver maximum iterations") = 50;
  
  
    system.add_variable("u",
  		      Utility::string_to_enum<Order>   (order),
  		      Utility::string_to_enum<FEFamily>(family));
  
    system.nonlinear_solver->residual = compute_residual;
    system.nonlinear_solver->jacobian = compute_jacobian;
  
    equation_systems.init();
  
    equation_systems.print_info();
  
    equation_systems.get_system("Laplace-Young").solve();
  
    std::cout << "Laplace-Young system solved at nonlinear iteration "
  	    << system.n_nonlinear_iterations()
  	    << " , final nonlinear residual norm: "
  	    << system.final_nonlinear_residual()
  	    << std::endl;
  
  #ifdef LIBMESH_HAVE_EXODUS_API
    ExodusII_IO (mesh).write_equation_systems ("out.e",
                                         equation_systems);
  #endif // #ifdef LIBMESH_HAVE_EXODUS_API
  #endif // #ifndef LIBMESH_ENABLE_AMR
  
    return 0;
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/miscellaneous/miscellaneous_ex3'
***************************************************************
* Running Example miscellaneous_ex3:
*  mpirun -np 4 example-devel -r 3 -o FIRST -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
 
Running /net/spark/workspace/roystgnr/libmesh/git/devel/examples/miscellaneous/miscellaneous_ex3/.libs/lt-example-devel -r 3 -o FIRST -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc

 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=225
    n_local_nodes()=66
  n_elem()=255
    n_local_elem()=70
    n_active_elem()=192
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "Laplace-Young"
    Type "NonlinearImplicit"
    Variables="u" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="FIRST", "THIRD" 
    n_dofs()=225
    n_local_dofs()=66
    n_constrained_dofs()=0
    n_local_constrained_dofs()=0
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 7.68444
      Average Off-Processor Bandwidth <= 0.853333
      Maximum  On-Processor Bandwidth <= 11
      Maximum Off-Processor Bandwidth <= 6
    DofMap Constraints
      Number of DoF Constraints = 0
      Number of Node Constraints = 0

  NL step  0, |residual|_2 = 2.000000e-01
  NL step  1, |residual|_2 = 4.433526e-03
  NL step  2, |residual|_2 = 2.163914e-04
  NL step  3, |residual|_2 = 1.157687e-05
  NL step  4, |residual|_2 = 6.567375e-07
  NL step  5, |residual|_2 = 3.849456e-08
  NL step  6, |residual|_2 = 2.293554e-09
Laplace-Young system solved at nonlinear iteration 6 , final nonlinear residual norm: 2.293554e-09

 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:50:53 2013                                                                          |
| OS:             Linux                                                                                             |
| HostName:       spark.ices.utexas.edu                                                                             |
| OS Release:     2.6.32-279.22.1.el6.x86_64                                                                        |
| OS Version:     #1 SMP Tue Feb 5 14:33:39 CST 2013                                                                |
| Machine:        x86_64                                                                                            |
| Username:       roystgnr                                                                                          |
| Configuration:  ../configure  '--enable-everything'                                                               |
|  'METHODS=devel'                                                                                                  |
|  '--prefix=/h2/roystgnr/libmesh-test'                                                                             |
|  'CXX=distcc /usr/bin/g++'                                                                                        |
|  'CC=distcc /usr/bin/gcc'                                                                                         |
|  'FC=distcc /usr/bin/gfortran'                                                                                    |
|  'F77=distcc /usr/bin/gfortran'                                                                                   |
|  'PETSC_DIR=/opt/apps/ossw/libraries/petsc/petsc-3.3-p2'                                                          |
|  'PETSC_ARCH=gcc-system-mkl-gf-10.3.12.361-mpich2-1.4.1p1-cxx-opt'                                                |
|  'SLEPC_DIR=/opt/apps/ossw/libraries/slepc/slepc-3.3-p2-petsc-3.3-p2-cxx-opt'                                     |
|  'TRILINOS_DIR=/opt/apps/ossw/libraries/trilinos/trilinos-10.12.2/sl6/gcc-system/mpich2-1.4.1p1/mkl-gf-10.3.12.361'|
|  'VTK_DIR=/opt/apps/ossw/libraries/vtk/vtk-5.10.0/sl6/gcc-system'                                                 |
|  'HDF5_DIR=/opt/apps/ossw/libraries/hdf5/hdf5-1.8.9/sl6/gcc-system'                                               |
 -------------------------------------------------------------------------------------------------------------------
 ----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=0.434914, Active time=0.237467                                                 |
 ----------------------------------------------------------------------------------------------------------------
| Event                              nCalls    Total Time  Avg Time    Total Time  Avg Time    % of Active Time  |
|                                              w/o Sub     w/o Sub     With Sub    With Sub    w/o S    With S   |
|----------------------------------------------------------------------------------------------------------------|
|                                                                                                                |
|                                                                                                                |
| DofMap                                                                                                         |
|   add_neighbors_to_send_list()     1         0.0003      0.000330    0.0006      0.000557    0.14     0.23     |
|   build_sparsity()                 1         0.0004      0.000359    0.0009      0.000946    0.15     0.40     |
|   create_dof_constraints()         1         0.0002      0.000229    0.0002      0.000229    0.10     0.10     |
|   distribute_dofs()                1         0.0007      0.000735    0.0020      0.002036    0.31     0.86     |
|   dof_indices()                    761       0.0044      0.000006    0.0044      0.000006    1.86     1.86     |
|   prepare_send_list()              1         0.0000      0.000006    0.0000      0.000006    0.00     0.00     |
|   reinit()                         1         0.0010      0.000952    0.0010      0.000952    0.40     0.40     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          1         0.0002      0.000236    0.0008      0.000777    0.10     0.33     |
|                                                                                                                |
| ExodusII_IO                                                                                                    |
|   write_nodal_data()               1         0.1547      0.154657    0.1547      0.154657    65.13    65.13    |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        722       0.0036      0.000005    0.0036      0.000005    1.51     1.51     |
|   init_shape_functions()           111       0.0003      0.000003    0.0003      0.000003    0.14     0.14     |
|   inverse_map()                    294       0.0011      0.000004    0.0011      0.000004    0.45     0.45     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             722       0.0025      0.000003    0.0025      0.000003    1.05     1.05     |
|   compute_face_map()               98        0.0008      0.000008    0.0019      0.000019    0.34     0.80     |
|   init_face_shape_functions()      7         0.0000      0.000004    0.0000      0.000004    0.01     0.01     |
|   init_reference_to_physical_map() 111       0.0008      0.000007    0.0008      0.000007    0.32     0.32     |
|                                                                                                                |
| LocationMap                                                                                                    |
|   find()                           756       0.0010      0.000001    0.0010      0.000001    0.40     0.40     |
|   init()                           3         0.0001      0.000023    0.0001      0.000023    0.03     0.03     |
|                                                                                                                |
| Mesh                                                                                                           |
|   find_neighbors()                 2         0.0011      0.000531    0.0012      0.000615    0.45     0.52     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   broadcast()                      1         0.0002      0.000241    0.0005      0.000548    0.10     0.23     |
|   compute_hilbert_indices()        3         0.0008      0.000277    0.0008      0.000277    0.35     0.35     |
|   find_global_indices()            3         0.0002      0.000083    0.0021      0.000702    0.10     0.89     |
|   parallel_sort()                  3         0.0004      0.000150    0.0006      0.000203    0.19     0.26     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         1         0.0001      0.000068    0.1556      0.155581    0.03     65.52    |
|                                                                                                                |
| MeshRefinement                                                                                                 |
|   _refine_elements()               3         0.0013      0.000445    0.0041      0.001358    0.56     1.72     |
|   add_point()                      756       0.0014      0.000002    0.0025      0.000003    0.61     1.04     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      2         0.0026      0.001279    0.0043      0.002138    1.08     1.80     |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      11        0.0001      0.000013    0.0002      0.000018    0.06     0.08     |
|   broadcast()                      27        0.0001      0.000005    0.0001      0.000003    0.05     0.03     |
|   max(bool)                        4         0.0001      0.000035    0.0001      0.000035    0.06     0.06     |
|   max(scalar)                      261       0.0013      0.000005    0.0013      0.000005    0.56     0.56     |
|   max(vector)                      62        0.0005      0.000008    0.0015      0.000024    0.21     0.62     |
|   min(bool)                        312       0.0014      0.000004    0.0014      0.000004    0.57     0.57     |
|   min(scalar)                      253       0.0019      0.000007    0.0019      0.000007    0.79     0.79     |
|   min(vector)                      62        0.0021      0.000034    0.0030      0.000048    0.89     1.26     |
|   probe()                          48        0.0001      0.000003    0.0001      0.000003    0.06     0.06     |
|   receive()                        48        0.0002      0.000003    0.0003      0.000006    0.06     0.13     |
|   send()                           48        0.0001      0.000002    0.0001      0.000002    0.04     0.04     |
|   send_receive()                   54        0.0003      0.000005    0.0008      0.000014    0.11     0.32     |
|   sum()                            36        0.0003      0.000007    0.0004      0.000011    0.11     0.17     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           48        0.0001      0.000001    0.0001      0.000001    0.03     0.03     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         2         0.0003      0.000168    0.0007      0.000329    0.14     0.28     |
|   set_parent_processor_ids()       2         0.0001      0.000061    0.0001      0.000061    0.05     0.05     |
|                                                                                                                |
| PetscNonlinearSolver                                                                                           |
|   jacobian()                       6         0.0096      0.001603    0.0166      0.002770    4.05     7.00     |
|   residual()                       7         0.0096      0.001366    0.0191      0.002734    4.03     8.06     |
|   solve()                          1         0.0290      0.028998    0.0648      0.064763    12.21    27.27    |
|                                                                                                                |
| System                                                                                                         |
|   solve()                          1         0.0000      0.000024    0.0648      0.064787    0.01     27.28    |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            5660      0.2375                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example miscellaneous_ex3:
*  mpirun -np 4 example-devel -r 3 -o FIRST -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
make[4]: Leaving directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/miscellaneous/miscellaneous_ex3'

Site Created By: libMesh Developers
Last modified: April 23 2013 04:19:30 UTC

Hosted By:
SourceForge.net Logo