The source file miscellaneous_ex5.C with comments:

Miscellaneous Example 5 - Interior Penalty Discontinuous Galerkin



By Lorenzo Botti

This example is based on Adaptivity Example 3, but uses an Interior Penalty Discontinuous Galerkin formulation.

        #include <iostream>
        
LibMesh include files.
        #include "libmesh/libmesh.h"
        #include "libmesh/mesh.h"
        #include "libmesh/equation_systems.h"
        #include "libmesh/mesh_data.h"
        #include "libmesh/mesh_generation.h"
        #include "libmesh/mesh_modification.h"
        #include "libmesh/elem.h"
        #include "libmesh/transient_system.h"
        #include "libmesh/fe.h"
        #include "libmesh/quadrature_gauss.h"
        #include "libmesh/dof_map.h"
        #include "libmesh/sparse_matrix.h"
        #include "libmesh/dense_matrix.h"
        #include "libmesh/dense_vector.h"
        #include "libmesh/dense_submatrix.h"
        #include "libmesh/dense_subvector.h"
        #include "libmesh/numeric_vector.h"
        #include "libmesh/linear_implicit_system.h"
        #include "libmesh/exodusII_io.h"
        #include "libmesh/fe_interface.h"
        #include "libmesh/getpot.h"
        #include "libmesh/mesh_refinement.h"
        #include "libmesh/error_vector.h"
        #include "libmesh/kelly_error_estimator.h"
        #include "libmesh/discontinuity_measure.h"
        #include "libmesh/string_to_enum.h"
        
        #include "libmesh/exact_solution.h"
#define QORDER TWENTYSIXTH

Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        Number exact_solution (const Point& p, const Parameters& parameters, const std::string&, const std::string&)
        {
          const Real x = p(0);
          const Real y = p(1);
          const Real z = p(2);
        
          if (parameters.get<bool>("singularity"))
          {
              Real theta = atan2(y,x);
        
              if (theta < 0)
                 theta += 2. * libMesh::pi;
        
              return pow(x*x + y*y, 1./3.)*sin(2./3.*theta) + z;
          }
          else
          {
              return cos(x) * exp(y) * (1. - z);
          }
        }
        
We now define the gradient of the exact solution, again being careful to obtain an angle from atan2 in the correct quadrant.

        Gradient exact_derivative(const Point& p,
                                  const Parameters& parameters,  // es parameters
                                  const std::string&,            // sys_name, not needed
                                  const std::string&)            // unk_name, not needed
        {
Gradient value to be returned.
          Gradient gradu;
        
x and y coordinates in space
          const Real x = p(0);
          const Real y = p(1);
          const Real z = p(2);
        
          if (parameters.get<bool>("singularity"))
          {
              libmesh_assert_not_equal_to (x, 0.);
        
For convenience...
              const Real tt = 2./3.;
              const Real ot = 1./3.;
        
The value of the radius, squared
              const Real r2 = x*x + y*y;
        
The boundary value, given by the exact solution, u_exact = r^(2/3)*sin(2*theta/3).
              Real theta = atan2(y,x);
        
Make sure 0 <= theta <= 2*pi
              if (theta < 0)
                theta += 2. * libMesh::pi;
        
du/dx
              gradu(0) = tt*x*pow(r2,-tt)*sin(tt*theta) - pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*y/x/x;
              gradu(1) = tt*y*pow(r2,-tt)*sin(tt*theta) + pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*1./x;
              gradu(2) = 1.;
          }
          else
          {
              gradu(0) = -sin(x) * exp(y) * (1. - z);
              gradu(1) = cos(x) * exp(y) * (1. - z);
              gradu(2) = -cos(x) * exp(y);
          }
          return gradu;
        }
        
We now define the matrix assembly function for the Laplace system. We need to first compute element volume matrices, and then take into account the boundary conditions and the flux integrals, which will be handled via an interior penalty method.
        void assemble_ellipticdg(EquationSystems& es, const std::string& system_name)
        {
          std::cout<<" assembling elliptic dg system... ";
          std::cout.flush();
        
It is a good idea to make sure we are assembling the proper system.
          libmesh_assert_equal_to (system_name, "EllipticDG");
        
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 & ellipticdg_system = es.get_system<LinearImplicitSystem> ("EllipticDG");
Get some parameters that we need during assembly
          const Real penalty = es.parameters.get<Real> ("penalty");
          std::string refinement_type = es.parameters.get<std::string> ("refinement");
        
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
          const DofMap & dof_map = ellipticdg_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 = ellipticdg_system.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));
          AutoPtr<FEBase> fe_elem_face(FEBase::build(dim, fe_type));
          AutoPtr<FEBase> fe_neighbor_face(FEBase::build(dim, fe_type));
        
Quadrature rules for numerical integration.
        #ifdef QORDER
          QGauss qrule (dim, QORDER);
        #else
          QGauss qrule (dim, fe_type.default_quadrature_order());
        #endif
          fe->attach_quadrature_rule (&qrule);
        
        #ifdef QORDER
          QGauss qface(dim-1, QORDER);
        #else
          QGauss qface(dim-1, fe_type.default_quadrature_order());
        #endif
        
Tell the finite element object to use our quadrature rule.
          fe_elem_face->attach_quadrature_rule(&qface);
          fe_neighbor_face->attach_quadrature_rule(&qface);
        
Here we define some references to cell-specific data that will be used to assemble the linear system. Data for interior volume integrals
          const std::vector<Real>& JxW = fe->get_JxW();
          const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
        
Data for surface integrals on the element boundary
          const std::vector<std::vector<Real> >&  phi_face = fe_elem_face->get_phi();
          const std::vector<std::vector<RealGradient> >& dphi_face = fe_elem_face->get_dphi();
          const std::vector<Real>& JxW_face = fe_elem_face->get_JxW();
          const std::vector<Point>& qface_normals = fe_elem_face->get_normals();
          const std::vector<Point>& qface_points = fe_elem_face->get_xyz();
        
Data for surface integrals on the neighbor boundary
          const std::vector<std::vector<Real> >&  phi_neighbor_face = fe_neighbor_face->get_phi();
          const std::vector<std::vector<RealGradient> >& dphi_neighbor_face = fe_neighbor_face->get_dphi();
        
Define data structures to contain the element interior 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;
        
Data structures to contain the element and neighbor boundary matrix contribution. This matrices will do the coupling beetwen the dofs of the element and those of his neighbors. Ken: matrix coupling elem and neighbor dofs
          DenseMatrix<Number> Kne;
          DenseMatrix<Number> Ken;
          DenseMatrix<Number> Kee;
          DenseMatrix<Number> Knn;
        
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 elements in the mesh. We will compute first the element interior matrix and right-hand-side contribution and then the element and neighbors boundary matrix contributions.
          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);
            const unsigned int n_dofs   = dof_indices.size();
        
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 matrix and right-hand side before summing them. We use the resize member here because the number of degrees of freedom might have changed from the last element.
            Ke.resize (n_dofs, n_dofs);
            Fe.resize (n_dofs);
        
Now we will build the element interior matrix. This involves a double loop to integrate the test funcions (i) against the trial functions (j).
            for (unsigned int qp=0; qp<qrule.n_points(); qp++)
            {
              for (unsigned int i=0; i<n_dofs; i++)
              {
                for (unsigned int j=0; j<n_dofs; j++)
                {
                   Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
                }
              }
            }
        
Now we adress boundary conditions. We consider Dirichlet bc imposed via the interior penalty method 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)
              {
Pointer to the element face
                fe_elem_face->reinit(elem, side);
        
                AutoPtr<Elem> elem_side (elem->build_side(side));
h elemet dimension to compute the interior penalty penalty parameter
                const unsigned int elem_b_order = static_cast<unsigned int> (fe_elem_face->get_order());
                const double h_elem = elem->volume()/elem_side->volume() * 1./pow(elem_b_order, 2.);
        
                for (unsigned int qp=0; qp<qface.n_points(); qp++)
                {
                  Number bc_value = exact_solution(qface_points[qp], es.parameters,"null","void");
                  for (unsigned int i=0; i<n_dofs; i++)
                  {
Matrix contribution
                     for (unsigned int j=0; j<n_dofs; j++)
                     {
                        Ke(i,j) += JxW_face[qp] * penalty/h_elem * phi_face[i][qp] * phi_face[j][qp]; // stability
                        Ke(i,j) -= JxW_face[qp] * (phi_face[i][qp] * (dphi_face[j][qp]*qface_normals[qp]) + phi_face[j][qp] * (dphi_face[i][qp]*qface_normals[qp])); // consistency
                     }
RHS contribution
                     Fe(i) += JxW_face[qp] * bc_value * penalty/h_elem * phi_face[i][qp];      // stability
                     Fe(i) -= JxW_face[qp] * dphi_face[i][qp] * (bc_value*qface_normals[qp]);     // consistency
                  }
                }
              }
If the element is not on a boundary of the domain we loop over his neighbors to compute the element and neighbor boundary matrix contributions
              else
              {
Store a pointer to the neighbor we are currently working on.
                const Elem* neighbor = elem->neighbor(side);
        
Get the global id of the element and the neighbor
                const unsigned int elem_id = elem->id();
                const unsigned int neighbor_id = neighbor->id();
        
If the neighbor has the same h level and is active perform integration only if our global id is bigger than our neighbor id. We don't want to compute twice the same contributions. If the neighbor has a different h level perform integration only if the neighbor is at a lower level.
                if ((neighbor->active() && (neighbor->level() == elem->level()) && (elem_id < neighbor_id)) || (neighbor->level() < elem->level()))
                {
Pointer to the element side
                  AutoPtr<Elem> elem_side (elem->build_side(side));
        
h dimension to compute the interior penalty penalty parameter
                  const unsigned int elem_b_order = static_cast<unsigned int>(fe_elem_face->get_order());
                  const unsigned int neighbor_b_order = static_cast<unsigned int>(fe_neighbor_face->get_order());
                  const double side_order = (elem_b_order + neighbor_b_order)/2.;
                  const double h_elem = (elem->volume()/elem_side->volume()) * 1./pow(side_order,2.);
        
The quadrature point locations on the neighbor side
                  std::vector<Point> qface_neighbor_point;
        
The quadrature point locations on the element side
                  std::vector<Point > qface_point;
        
Reinitialize shape functions on the element side
                  fe_elem_face->reinit(elem, side);
        
Get the physical locations of the element quadrature points
                  qface_point = fe_elem_face->get_xyz();
        
Find their locations on the neighbor
                  unsigned int side_neighbor = neighbor->which_neighbor_am_i(elem);
                  if (refinement_type == "p")
                    fe_neighbor_face->side_map (neighbor, elem_side.get(), side_neighbor, qface.get_points(), qface_neighbor_point);
        	  else
                    FEInterface::inverse_map (elem->dim(), fe->get_fe_type(), neighbor, qface_point, qface_neighbor_point);
Calculate the neighbor element shape functions at those locations
                  fe_neighbor_face->reinit(neighbor, &qface_neighbor_point);
        
Get the degree of freedom indices for the neighbor. These define where in the global matrix this neighbor will contribute to.
                  std::vector<dof_id_type> neighbor_dof_indices;
                  dof_map.dof_indices (neighbor, neighbor_dof_indices);
                  const unsigned int n_neighbor_dofs = neighbor_dof_indices.size();
        
Zero the element and neighbor side matrix before summing them. We use the resize member here because the number of degrees of freedom might have changed from the last element or neighbor. Note that Kne and Ken are not square matrices if neighbor and element have a different p level
                  Kne.resize (n_neighbor_dofs, n_dofs);
                  Ken.resize (n_dofs, n_neighbor_dofs);
                  Kee.resize (n_dofs, n_dofs);
                  Knn.resize (n_neighbor_dofs, n_neighbor_dofs);
        
Now we will build the element and neighbor boundary matrices. This involves a double loop to integrate the test funcions (i) against the trial functions (j).
                  for (unsigned int qp=0; qp<qface.n_points(); qp++)
                  {
Kee Matrix. Integrate the element test function i against the element test function j
                    for (unsigned int i=0; i<n_dofs; i++)
                    {
                      for (unsigned int j=0; j<n_dofs; j++)
                      {
                         Kee(i,j) -= 0.5 * JxW_face[qp] * (phi_face[j][qp]*(qface_normals[qp]*dphi_face[i][qp]) + phi_face[i][qp]*(qface_normals[qp]*dphi_face[j][qp])); // consistency
                         Kee(i,j) += JxW_face[qp] * penalty/h_elem * phi_face[j][qp]*phi_face[i][qp];  // stability
                      }
                    }
        
Knn Matrix. Integrate the neighbor test function i against the neighbor test function j
                    for (unsigned int i=0; i<n_neighbor_dofs; i++)
                    {
                      for (unsigned int j=0; j<n_neighbor_dofs; j++)
                      {
                         Knn(i,j) += 0.5 * JxW_face[qp] * (phi_neighbor_face[j][qp]*(qface_normals[qp]*dphi_neighbor_face[i][qp]) + phi_neighbor_face[i][qp]*(qface_normals[qp]*dphi_neighbor_face[j][qp])); // consistency
                         Knn(i,j) += JxW_face[qp] * penalty/h_elem * phi_neighbor_face[j][qp]*phi_neighbor_face[i][qp]; // stability
                      }
                    }
        
Kne Matrix. Integrate the neighbor test function i against the element test function j
                    for (unsigned int i=0; i<n_neighbor_dofs; i++)
                    {
                      for (unsigned int j=0; j<n_dofs; j++)
                      {
                         Kne(i,j) += 0.5 * JxW_face[qp] * (phi_neighbor_face[i][qp]*(qface_normals[qp]*dphi_face[j][qp]) - phi_face[j][qp]*(qface_normals[qp]*dphi_neighbor_face[i][qp])); // consistency
                         Kne(i,j) -= JxW_face[qp] * penalty/h_elem * phi_face[j][qp]*phi_neighbor_face[i][qp];  // stability
                      }
                    }
        
Ken Matrix. Integrate the element test function i against the neighbor test function j
                    for (unsigned int i=0; i<n_dofs; i++)
                    {
                      for (unsigned int j=0; j<n_neighbor_dofs; j++)
                      {
                         Ken(i,j) += 0.5 * JxW_face[qp] * (phi_neighbor_face[j][qp]*(qface_normals[qp]*dphi_face[i][qp]) - phi_face[i][qp]*(qface_normals[qp]*dphi_neighbor_face[j][qp]));  // consistency
                         Ken(i,j) -= JxW_face[qp] * penalty/h_elem * phi_face[i][qp]*phi_neighbor_face[j][qp];  // stability
                      }
                    }
                  }
        
The element and neighbor boundary matrix are now built for this side. Add them to the global matrix The \p SparseMatrix::add_matrix() members do this for us.
                  ellipticdg_system.matrix->add_matrix(Kne,neighbor_dof_indices,dof_indices);
                  ellipticdg_system.matrix->add_matrix(Ken,dof_indices,neighbor_dof_indices);
                  ellipticdg_system.matrix->add_matrix(Kee,dof_indices);
                  ellipticdg_system.matrix->add_matrix(Knn,neighbor_dof_indices);
                }
              }
            }
The element interior matrix and right-hand-side are now built for this element. Add them to the global matrix and right-hand-side vector. The \p SparseMatrix::add_matrix() and \p NumericVector::add_vector() members do this for us.
            ellipticdg_system.matrix->add_matrix(Ke, dof_indices);
            ellipticdg_system.rhs->add_vector(Fe, dof_indices);
          }
        
          std::cout << "done" << std::endl;
        }
        
        
        
        int main (int argc, char** argv)
        {
          LibMeshInit init(argc, argv);
        
Parse the input file
          GetPot input_file("miscellaneous_ex5.in");
        
Read in parameters from the input file
          const unsigned int adaptive_refinement_steps = input_file("max_adaptive_r_steps",3);
          const unsigned int uniform_refinement_steps  = input_file("uniform_h_r_steps",3);
          const Real refine_fraction                   = input_file("refine_fraction",0.5);
          const Real coarsen_fraction                  = input_file("coarsen_fraction",0.);
          const unsigned int max_h_level               = input_file("max_h_level", 10);
          const std::string refinement_type            = input_file("refinement_type","p");
          Order p_order                                = static_cast<Order>(input_file("p_order",1));
          const std::string element_type               = input_file("element_type", "tensor");
          const Real penalty                           = input_file("ip_penalty", 10.);
          const bool singularity                       = input_file("singularity", true);
          const unsigned int dim                       = input_file("dimension", 3);
        
Skip higher-dimensional examples on a lower-dimensional libMesh build
          libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");
        
Skip adaptive examples on a non-adaptive libMesh build
        #ifndef LIBMESH_ENABLE_AMR
          libmesh_example_assert(false, "--enable-amr");
        #else
        
Create a mesh, with dimension to be overridden later, distributed across the default MPI communicator.
          Mesh mesh(init.comm());
        
          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);
        
Mesh Refinement object
          MeshRefinement mesh_refinement(mesh);
          mesh_refinement.refine_fraction() = refine_fraction;
          mesh_refinement.coarsen_fraction() = coarsen_fraction;
          mesh_refinement.max_h_level() = max_h_level;
Do uniform refinement
          for (unsigned int rstep=0; rstep<uniform_refinement_steps; rstep++)
          {
             mesh_refinement.uniformly_refine(1);
          }
        
Crate an equation system object
          EquationSystems equation_system (mesh);
        
Set parameters for the equation system and the solver
          equation_system.parameters.set<Real>("linear solver tolerance") = TOLERANCE * TOLERANCE;
          equation_system.parameters.set<unsigned int>("linear solver maximum iterations") = 1000;
          equation_system.parameters.set<Real>("penalty") = penalty;
          equation_system.parameters.set<bool>("singularity") = singularity;
          equation_system.parameters.set<std::string>("refinement") = refinement_type;
        
Create a system named ellipticdg
          LinearImplicitSystem& ellipticdg_system = equation_system.add_system<LinearImplicitSystem> ("EllipticDG");
        
Add a variable "u" to "ellipticdg" using the p_order specified in the config file
          if( on_command_line( "element_type" ) )
            {
              std::string fe_str = command_line_value( std::string("element_type"),
        					       std::string("MONOMIAL") );
              if( fe_str != "MONOMIAL" || fe_str != "XYZ" )
        	{
        	  std::cerr << "Error: This example must be run with MONOMIAL or XYZ element types." << std::endl;
        	  libmesh_error();
        	}
              ellipticdg_system.add_variable ("u", p_order, Utility::string_to_enum<FEFamily>(fe_str) );
             }
          else
            ellipticdg_system.add_variable ("u", p_order, MONOMIAL);
        
Give the system a pointer to the matrix assembly function
          ellipticdg_system.attach_assemble_function (assemble_ellipticdg);
        
Initialize the data structures for the equation system
          equation_system.init();
        
Construct ExactSolution object and attach solution functions
          ExactSolution exact_sol(equation_system);
          exact_sol.attach_exact_value(exact_solution);
          exact_sol.attach_exact_deriv(exact_derivative);
        
A refinement loop.
          for (unsigned int rstep=0; rstep<adaptive_refinement_steps; ++rstep)
          {
            std::cout << "  Beginning Solve " << rstep << std::endl;
            std::cout << "Number of elements: " << mesh.n_elem() << std::endl;
        
Solve the system
            ellipticdg_system.solve();
        
            std::cout << "System has: " << equation_system.n_active_dofs()
                      << " degrees of freedom."
                      << std::endl;
        
            std::cout << "Linear solver converged at step: "
                      << ellipticdg_system.n_linear_iterations()
                      << ", final residual: "
                      << ellipticdg_system.final_linear_residual()
                      << std::endl;
        
Compute the error
            exact_sol.compute_error("EllipticDG", "u");
        
Print out the error values
            std::cout << "L2-Error is: "
                      << exact_sol.l2_error("EllipticDG", "u")
                      << std::endl;
        
Possibly refine the mesh
            if (rstep+1 < adaptive_refinement_steps)
            {
The ErrorVector is a particular StatisticsVector for computing error information on a finite element mesh.
              ErrorVector error;
        
The discontinuity error estimator evaluate the jump of the solution on elements faces
              DiscontinuityMeasure error_estimator;
              error_estimator.estimate_error(ellipticdg_system,error);
        
Take the error in error and decide which elements will be coarsened or refined
              mesh_refinement.flag_elements_by_error_fraction(error);
              if (refinement_type == "p")
                mesh_refinement.switch_h_to_p_refinement();
              if (refinement_type == "hp")
                mesh_refinement.add_p_to_h_refinement();
Refine and coarsen the flagged elements
              mesh_refinement.refine_and_coarsen_elements();
              equation_system.reinit();
            }
          }
        
Write out the solution After solving the system write the solution to a ExodusII-formatted plot file.
        #ifdef LIBMESH_HAVE_EXODUS_API
          ExodusII_IO (mesh).write_discontinuous_exodusII("lshaped_dg.e",equation_system);
        #endif
        
        #endif // #ifndef LIBMESH_ENABLE_AMR
        
All done.
          return 0;
        }



The source file miscellaneous_ex5.C without comments:

 
  
  #include <iostream>
  
  #include "libmesh/libmesh.h"
  #include "libmesh/mesh.h"
  #include "libmesh/equation_systems.h"
  #include "libmesh/mesh_data.h"
  #include "libmesh/mesh_generation.h"
  #include "libmesh/mesh_modification.h"
  #include "libmesh/elem.h"
  #include "libmesh/transient_system.h"
  #include "libmesh/fe.h"
  #include "libmesh/quadrature_gauss.h"
  #include "libmesh/dof_map.h"
  #include "libmesh/sparse_matrix.h"
  #include "libmesh/dense_matrix.h"
  #include "libmesh/dense_vector.h"
  #include "libmesh/dense_submatrix.h"
  #include "libmesh/dense_subvector.h"
  #include "libmesh/numeric_vector.h"
  #include "libmesh/linear_implicit_system.h"
  #include "libmesh/exodusII_io.h"
  #include "libmesh/fe_interface.h"
  #include "libmesh/getpot.h"
  #include "libmesh/mesh_refinement.h"
  #include "libmesh/error_vector.h"
  #include "libmesh/kelly_error_estimator.h"
  #include "libmesh/discontinuity_measure.h"
  #include "libmesh/string_to_enum.h"
  
  #include "libmesh/exact_solution.h"
  
  using namespace libMesh;
  
  Number exact_solution (const Point& p, const Parameters& parameters, const std::string&, const std::string&)
  {
    const Real x = p(0);
    const Real y = p(1);
    const Real z = p(2);
  
    if (parameters.get<bool>("singularity"))
    {
        Real theta = atan2(y,x);
  
        if (theta < 0)
           theta += 2. * libMesh::pi;
  
        return pow(x*x + y*y, 1./3.)*sin(2./3.*theta) + z;
    }
    else
    {
        return cos(x) * exp(y) * (1. - z);
    }
  }
  
  
  Gradient exact_derivative(const Point& p,
                            const Parameters& parameters,  // es parameters
                            const std::string&,            // sys_name, not needed
                            const std::string&)            // unk_name, not needed
  {
    Gradient gradu;
  
    const Real x = p(0);
    const Real y = p(1);
    const Real z = p(2);
  
    if (parameters.get<bool>("singularity"))
    {
        libmesh_assert_not_equal_to (x, 0.);
  
        const Real tt = 2./3.;
        const Real ot = 1./3.;
  
        const Real r2 = x*x + y*y;
  
        Real theta = atan2(y,x);
  
        if (theta < 0)
          theta += 2. * libMesh::pi;
  
        gradu(0) = tt*x*pow(r2,-tt)*sin(tt*theta) - pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*y/x/x;
        gradu(1) = tt*y*pow(r2,-tt)*sin(tt*theta) + pow(r2,ot)*cos(tt*theta)*tt/(1.+y*y/x/x)*1./x;
        gradu(2) = 1.;
    }
    else
    {
        gradu(0) = -sin(x) * exp(y) * (1. - z);
        gradu(1) = cos(x) * exp(y) * (1. - z);
        gradu(2) = -cos(x) * exp(y);
    }
    return gradu;
  }
  
  void assemble_ellipticdg(EquationSystems& es, const std::string& system_name)
  {
    std::cout<<" assembling elliptic dg system... ";
    std::cout.flush();
  
    libmesh_assert_equal_to (system_name, "EllipticDG");
  
    const MeshBase& mesh = es.get_mesh();
    const unsigned int dim = mesh.mesh_dimension();
  
    LinearImplicitSystem & ellipticdg_system = es.get_system<LinearImplicitSystem> ("EllipticDG");
    const Real penalty = es.parameters.get<Real> ("penalty");
    std::string refinement_type = es.parameters.get<std::string> ("refinement");
  
    const DofMap & dof_map = ellipticdg_system.get_dof_map();
  
    FEType fe_type = ellipticdg_system.variable_type(0);
  
    AutoPtr<FEBase> fe  (FEBase::build(dim, fe_type));
    AutoPtr<FEBase> fe_elem_face(FEBase::build(dim, fe_type));
    AutoPtr<FEBase> fe_neighbor_face(FEBase::build(dim, fe_type));
  
  #ifdef QORDER
    QGauss qrule (dim, QORDER);
  #else
    QGauss qrule (dim, fe_type.default_quadrature_order());
  #endif
    fe->attach_quadrature_rule (&qrule);
  
  #ifdef QORDER
    QGauss qface(dim-1, QORDER);
  #else
    QGauss qface(dim-1, fe_type.default_quadrature_order());
  #endif
  
    fe_elem_face->attach_quadrature_rule(&qface);
    fe_neighbor_face->attach_quadrature_rule(&qface);
  
    const std::vector<Real>& JxW = fe->get_JxW();
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
    const std::vector<std::vector<Real> >&  phi_face = fe_elem_face->get_phi();
    const std::vector<std::vector<RealGradient> >& dphi_face = fe_elem_face->get_dphi();
    const std::vector<Real>& JxW_face = fe_elem_face->get_JxW();
    const std::vector<Point>& qface_normals = fe_elem_face->get_normals();
    const std::vector<Point>& qface_points = fe_elem_face->get_xyz();
  
    const std::vector<std::vector<Real> >&  phi_neighbor_face = fe_neighbor_face->get_phi();
    const std::vector<std::vector<RealGradient> >& dphi_neighbor_face = fe_neighbor_face->get_dphi();
  
    DenseMatrix<Number> Ke;
    DenseVector<Number> Fe;
  
    DenseMatrix<Number> Kne;
    DenseMatrix<Number> Ken;
    DenseMatrix<Number> Kee;
    DenseMatrix<Number> Knn;
  
    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);
      const unsigned int n_dofs   = dof_indices.size();
  
      fe->reinit  (elem);
  
      Ke.resize (n_dofs, n_dofs);
      Fe.resize (n_dofs);
  
      for (unsigned int qp=0; qp<qrule.n_points(); qp++)
      {
        for (unsigned int i=0; i<n_dofs; i++)
        {
          for (unsigned int j=0; j<n_dofs; j++)
          {
             Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
          }
        }
      }
  
      for (unsigned int side=0; side<elem->n_sides(); side++)
      {
        if (elem->neighbor(side) == NULL)
        {
          fe_elem_face->reinit(elem, side);
  
          AutoPtr<Elem> elem_side (elem->build_side(side));
          const unsigned int elem_b_order = static_cast<unsigned int> (fe_elem_face->get_order());
          const double h_elem = elem->volume()/elem_side->volume() * 1./pow(elem_b_order, 2.);
  
          for (unsigned int qp=0; qp<qface.n_points(); qp++)
          {
            Number bc_value = exact_solution(qface_points[qp], es.parameters,"null","void");
            for (unsigned int i=0; i<n_dofs; i++)
            {
               for (unsigned int j=0; j<n_dofs; j++)
               {
                  Ke(i,j) += JxW_face[qp] * penalty/h_elem * phi_face[i][qp] * phi_face[j][qp]; // stability
                  Ke(i,j) -= JxW_face[qp] * (phi_face[i][qp] * (dphi_face[j][qp]*qface_normals[qp]) + phi_face[j][qp] * (dphi_face[i][qp]*qface_normals[qp])); // consistency
               }
               Fe(i) += JxW_face[qp] * bc_value * penalty/h_elem * phi_face[i][qp];      // stability
               Fe(i) -= JxW_face[qp] * dphi_face[i][qp] * (bc_value*qface_normals[qp]);     // consistency
            }
          }
        }
        else
        {
          const Elem* neighbor = elem->neighbor(side);
  
          const unsigned int elem_id = elem->id();
          const unsigned int neighbor_id = neighbor->id();
  
          if ((neighbor->active() && (neighbor->level() == elem->level()) && (elem_id < neighbor_id)) || (neighbor->level() < elem->level()))
          {
            AutoPtr<Elem> elem_side (elem->build_side(side));
  
            const unsigned int elem_b_order = static_cast<unsigned int>(fe_elem_face->get_order());
            const unsigned int neighbor_b_order = static_cast<unsigned int>(fe_neighbor_face->get_order());
            const double side_order = (elem_b_order + neighbor_b_order)/2.;
            const double h_elem = (elem->volume()/elem_side->volume()) * 1./pow(side_order,2.);
  
            std::vector<Point> qface_neighbor_point;
  
            std::vector<Point > qface_point;
  
            fe_elem_face->reinit(elem, side);
  
            qface_point = fe_elem_face->get_xyz();
  
  	  unsigned int side_neighbor = neighbor->which_neighbor_am_i(elem);
            if (refinement_type == "p")
              fe_neighbor_face->side_map (neighbor, elem_side.get(), side_neighbor, qface.get_points(), qface_neighbor_point);
  	  else
              FEInterface::inverse_map (elem->dim(), fe->get_fe_type(), neighbor, qface_point, qface_neighbor_point);
            fe_neighbor_face->reinit(neighbor, &qface_neighbor_point);
  
            std::vector<dof_id_type> neighbor_dof_indices;
            dof_map.dof_indices (neighbor, neighbor_dof_indices);
            const unsigned int n_neighbor_dofs = neighbor_dof_indices.size();
  
            Kne.resize (n_neighbor_dofs, n_dofs);
            Ken.resize (n_dofs, n_neighbor_dofs);
            Kee.resize (n_dofs, n_dofs);
            Knn.resize (n_neighbor_dofs, n_neighbor_dofs);
  
            for (unsigned int qp=0; qp<qface.n_points(); qp++)
            {
              for (unsigned int i=0; i<n_dofs; i++)
              {
                for (unsigned int j=0; j<n_dofs; j++)
                {
                   Kee(i,j) -= 0.5 * JxW_face[qp] * (phi_face[j][qp]*(qface_normals[qp]*dphi_face[i][qp]) + phi_face[i][qp]*(qface_normals[qp]*dphi_face[j][qp])); // consistency
                   Kee(i,j) += JxW_face[qp] * penalty/h_elem * phi_face[j][qp]*phi_face[i][qp];  // stability
                }
              }
  
              for (unsigned int i=0; i<n_neighbor_dofs; i++)
              {
                for (unsigned int j=0; j<n_neighbor_dofs; j++)
                {
                   Knn(i,j) += 0.5 * JxW_face[qp] * (phi_neighbor_face[j][qp]*(qface_normals[qp]*dphi_neighbor_face[i][qp]) + phi_neighbor_face[i][qp]*(qface_normals[qp]*dphi_neighbor_face[j][qp])); // consistency
                   Knn(i,j) += JxW_face[qp] * penalty/h_elem * phi_neighbor_face[j][qp]*phi_neighbor_face[i][qp]; // stability
                }
              }
  
              for (unsigned int i=0; i<n_neighbor_dofs; i++)
              {
                for (unsigned int j=0; j<n_dofs; j++)
                {
                   Kne(i,j) += 0.5 * JxW_face[qp] * (phi_neighbor_face[i][qp]*(qface_normals[qp]*dphi_face[j][qp]) - phi_face[j][qp]*(qface_normals[qp]*dphi_neighbor_face[i][qp])); // consistency
                   Kne(i,j) -= JxW_face[qp] * penalty/h_elem * phi_face[j][qp]*phi_neighbor_face[i][qp];  // stability
                }
              }
  
              for (unsigned int i=0; i<n_dofs; i++)
              {
                for (unsigned int j=0; j<n_neighbor_dofs; j++)
                {
                   Ken(i,j) += 0.5 * JxW_face[qp] * (phi_neighbor_face[j][qp]*(qface_normals[qp]*dphi_face[i][qp]) - phi_face[i][qp]*(qface_normals[qp]*dphi_neighbor_face[j][qp]));  // consistency
                   Ken(i,j) -= JxW_face[qp] * penalty/h_elem * phi_face[i][qp]*phi_neighbor_face[j][qp];  // stability
                }
              }
            }
  
            ellipticdg_system.matrix->add_matrix(Kne,neighbor_dof_indices,dof_indices);
            ellipticdg_system.matrix->add_matrix(Ken,dof_indices,neighbor_dof_indices);
            ellipticdg_system.matrix->add_matrix(Kee,dof_indices);
            ellipticdg_system.matrix->add_matrix(Knn,neighbor_dof_indices);
          }
        }
      }
      ellipticdg_system.matrix->add_matrix(Ke, dof_indices);
      ellipticdg_system.rhs->add_vector(Fe, dof_indices);
    }
  
    std::cout << "done" << std::endl;
  }
  
  
  
  int main (int argc, char** argv)
  {
    LibMeshInit init(argc, argv);
  
    GetPot input_file("miscellaneous_ex5.in");
  
    const unsigned int adaptive_refinement_steps = input_file("max_adaptive_r_steps",3);
    const unsigned int uniform_refinement_steps  = input_file("uniform_h_r_steps",3);
    const Real refine_fraction                   = input_file("refine_fraction",0.5);
    const Real coarsen_fraction                  = input_file("coarsen_fraction",0.);
    const unsigned int max_h_level               = input_file("max_h_level", 10);
    const std::string refinement_type            = input_file("refinement_type","p");
    Order p_order                                = static_cast<Order>(input_file("p_order",1));
    const std::string element_type               = input_file("element_type", "tensor");
    const Real penalty                           = input_file("ip_penalty", 10.);
    const bool singularity                       = input_file("singularity", true);
    const unsigned int dim                       = input_file("dimension", 3);
  
    libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");
  
  #ifndef LIBMESH_ENABLE_AMR
    libmesh_example_assert(false, "--enable-amr");
  #else
  
    Mesh mesh(init.comm());
  
    if (dim == 1)
      MeshTools::Generation::build_line(mesh,1,-1.,0.);
    else if (dim == 2)
      mesh.read("lshaped.xda");
    else
      mesh.read ("lshaped3D.xda");
  
    if (element_type == "simplex")
      MeshTools::Modification::all_tri(mesh);
  
    MeshRefinement mesh_refinement(mesh);
    mesh_refinement.refine_fraction() = refine_fraction;
    mesh_refinement.coarsen_fraction() = coarsen_fraction;
    mesh_refinement.max_h_level() = max_h_level;
    for (unsigned int rstep=0; rstep<uniform_refinement_steps; rstep++)
    {
       mesh_refinement.uniformly_refine(1);
    }
  
    EquationSystems equation_system (mesh);
  
    equation_system.parameters.set<Real>("linear solver tolerance") = TOLERANCE * TOLERANCE;
    equation_system.parameters.set<unsigned int>("linear solver maximum iterations") = 1000;
    equation_system.parameters.set<Real>("penalty") = penalty;
    equation_system.parameters.set<bool>("singularity") = singularity;
    equation_system.parameters.set<std::string>("refinement") = refinement_type;
  
    LinearImplicitSystem& ellipticdg_system = equation_system.add_system<LinearImplicitSystem> ("EllipticDG");
  
    if( on_command_line( "element_type" ) )
      {
        std::string fe_str = command_line_value( std::string("element_type"),
  					       std::string("MONOMIAL") );
        if( fe_str != "MONOMIAL" || fe_str != "XYZ" )
  	{
  	  std::cerr << "Error: This example must be run with MONOMIAL or XYZ element types." << std::endl;
  	  libmesh_error();
  	}
        ellipticdg_system.add_variable ("u", p_order, Utility::string_to_enum<FEFamily>(fe_str) );
       }
    else
      ellipticdg_system.add_variable ("u", p_order, MONOMIAL);
  
    ellipticdg_system.attach_assemble_function (assemble_ellipticdg);
  
    equation_system.init();
  
    ExactSolution exact_sol(equation_system);
    exact_sol.attach_exact_value(exact_solution);
    exact_sol.attach_exact_deriv(exact_derivative);
  
    for (unsigned int rstep=0; rstep<adaptive_refinement_steps; ++rstep)
    {
      std::cout << "  Beginning Solve " << rstep << std::endl;
      std::cout << "Number of elements: " << mesh.n_elem() << std::endl;
  
      ellipticdg_system.solve();
  
      std::cout << "System has: " << equation_system.n_active_dofs()
                << " degrees of freedom."
                << std::endl;
  
      std::cout << "Linear solver converged at step: "
                << ellipticdg_system.n_linear_iterations()
                << ", final residual: "
                << ellipticdg_system.final_linear_residual()
                << std::endl;
  
      exact_sol.compute_error("EllipticDG", "u");
  
      std::cout << "L2-Error is: "
                << exact_sol.l2_error("EllipticDG", "u")
                << std::endl;
  
      if (rstep+1 < adaptive_refinement_steps)
      {
        ErrorVector error;
  
        DiscontinuityMeasure error_estimator;
        error_estimator.estimate_error(ellipticdg_system,error);
  
        mesh_refinement.flag_elements_by_error_fraction(error);
        if (refinement_type == "p")
          mesh_refinement.switch_h_to_p_refinement();
        if (refinement_type == "hp")
          mesh_refinement.add_p_to_h_refinement();
        mesh_refinement.refine_and_coarsen_elements();
        equation_system.reinit();
      }
    }
  
  #ifdef LIBMESH_HAVE_EXODUS_API
    ExodusII_IO (mesh).write_discontinuous_exodusII("lshaped_dg.e",equation_system);
  #endif
  
  #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_ex5'
***************************************************************
* Running Example miscellaneous_ex5:
*  mpirun -np 4 example-devel  -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
 
  Beginning Solve 0
Number of elements: 219
 assembling elliptic dg system... done
System has: 768 degrees of freedom.
Linear solver converged at step: 30, final residual: 9.16568e-12
L2-Error is: 0.00666744
  Beginning Solve 1
Number of elements: 827
 assembling elliptic dg system... done
System has: 2896 degrees of freedom.
Linear solver converged at step: 41, final residual: 1.2975e-11
L2-Error is: 0.00264921
  Beginning Solve 2
Number of elements: 3003
 assembling elliptic dg system... done
System has: 10512 degrees of freedom.
Linear solver converged at step: 50, final residual: 3.07715e-11
L2-Error is: 0.0016323

 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:51:17 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=2.51885, Active time=2.35228                                                        |
 ---------------------------------------------------------------------------------------------------------------------
| 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()          3         0.0054      0.001790    0.0082      0.002743    0.23     0.35     |
|   build_sparsity()                      3         0.0515      0.017158    0.1783      0.059427    2.19     7.58     |
|   create_dof_constraints()              3         0.0183      0.006114    0.0392      0.013060    0.78     1.67     |
|   distribute_dofs()                     3         0.0119      0.003981    0.0376      0.012527    0.51     1.60     |
|   dof_indices()                         30553     0.1827      0.000006    0.1827      0.000006    7.77     7.77     |
|   old_dof_indices()                     1776      0.0103      0.000006    0.0103      0.000006    0.44     0.44     |
|   prepare_send_list()                   3         0.0001      0.000032    0.0001      0.000032    0.00     0.00     |
|   reinit()                              3         0.0187      0.006222    0.0187      0.006222    0.79     0.79     |
|                                                                                                                     |
| EquationSystems                                                                                                     |
|   build_discontinuous_solution_vector() 1         0.0119      0.011941    0.0277      0.027704    0.51     1.18     |
|                                                                                                                     |
| ExodusII_IO                                                                                                         |
|   write_nodal_data_discontinuous()      1         0.0958      0.095761    0.0958      0.095761    4.07     4.07     |
|                                                                                                                     |
| FE                                                                                                                  |
|   compute_shape_functions()             10017     0.0397      0.000004    0.0397      0.000004    1.69     1.69     |
|   init_shape_functions()                8247      0.0281      0.000003    0.0281      0.000003    1.19     1.19     |
|   inverse_map()                         24404     0.1824      0.000007    0.1824      0.000007    7.75     7.75     |
|                                                                                                                     |
| FEMap                                                                                                               |
|   compute_affine_map()                  10017     0.0542      0.000005    0.0542      0.000005    2.31     2.31     |
|   compute_face_map()                    3612      0.0184      0.000005    0.0184      0.000005    0.78     0.78     |
|   init_face_shape_functions()           5         0.0001      0.000013    0.0001      0.000013    0.00     0.00     |
|   init_reference_to_physical_map()      8247      0.1352      0.000016    0.1352      0.000016    5.75     5.75     |
|                                                                                                                     |
| JumpErrorEstimator                                                                                                  |
|   estimate_error()                      2         0.0182      0.009118    0.0905      0.045270    0.78     3.85     |
|                                                                                                                     |
| LocationMap                                                                                                         |
|   find()                                22344     0.0257      0.000001    0.0257      0.000001    1.09     1.09     |
|   init()                                6         0.0028      0.000464    0.0028      0.000464    0.12     0.12     |
|                                                                                                                     |
| Mesh                                                                                                                |
|   contract()                            2         0.0007      0.000342    0.0017      0.000828    0.03     0.07     |
|   find_neighbors()                      5         0.0257      0.005146    0.0315      0.006307    1.09     1.34     |
|   renumber_nodes_and_elem()             12        0.0032      0.000271    0.0032      0.000271    0.14     0.14     |
|                                                                                                                     |
| MeshCommunication                                                                                                   |
|   compute_hilbert_indices()             5         0.0157      0.003131    0.0157      0.003131    0.67     0.67     |
|   find_global_indices()                 5         0.0027      0.000535    0.0298      0.005960    0.11     1.27     |
|   parallel_sort()                       5         0.0023      0.000465    0.0055      0.001092    0.10     0.23     |
|                                                                                                                     |
| MeshRefinement                                                                                                      |
|   _coarsen_elements()                   4         0.0010      0.000256    0.0010      0.000262    0.04     0.04     |
|   _refine_elements()                    6         0.0381      0.006347    0.1178      0.019632    1.62     5.01     |
|   add_point()                           22344     0.0411      0.000002    0.0694      0.000003    1.75     2.95     |
|   make_coarsening_compatible()          8         0.0159      0.001992    0.0159      0.001992    0.68     0.68     |
|   make_flags_parallel_consistent()      6         0.0065      0.001081    0.0084      0.001397    0.28     0.36     |
|   make_refinement_compatible()          8         0.0009      0.000115    0.0010      0.000119    0.04     0.04     |
|                                                                                                                     |
| MetisPartitioner                                                                                                    |
|   partition()                           5         0.0916      0.018315    0.1281      0.025622    3.89     5.45     |
|                                                                                                                     |
| Parallel                                                                                                            |
|   allgather()                           19        0.0030      0.000156    0.0036      0.000190    0.13     0.15     |
|   broadcast()                           20        0.0001      0.000007    0.0001      0.000004    0.01     0.00     |
|   max(bool)                             23        0.0078      0.000337    0.0078      0.000337    0.33     0.33     |
|   max(scalar)                           558       0.0151      0.000027    0.0151      0.000027    0.64     0.64     |
|   max(vector)                           135       0.0035      0.000026    0.0129      0.000095    0.15     0.55     |
|   min(bool)                             697       0.0152      0.000022    0.0152      0.000022    0.64     0.64     |
|   min(scalar)                           547       0.0376      0.000069    0.0376      0.000069    1.60     1.60     |
|   min(vector)                           135       0.0037      0.000027    0.0134      0.000099    0.16     0.57     |
|   probe()                               192       0.0081      0.000042    0.0081      0.000042    0.34     0.34     |
|   receive()                             192       0.0007      0.000004    0.0088      0.000046    0.03     0.37     |
|   send()                                192       0.0004      0.000002    0.0004      0.000002    0.02     0.02     |
|   send_receive()                        202       0.0009      0.000005    0.0104      0.000052    0.04     0.44     |
|   sum()                                 25        0.0023      0.000094    0.0146      0.000583    0.10     0.62     |
|                                                                                                                     |
| Parallel::Request                                                                                                   |
|   wait()                                192       0.0002      0.000001    0.0002      0.000001    0.01     0.01     |
|                                                                                                                     |
| Partitioner                                                                                                         |
|   set_node_processor_ids()              6         0.0055      0.000924    0.0166      0.002773    0.24     0.71     |
|   set_parent_processor_ids()            5         0.0027      0.000536    0.0027      0.000536    0.11     0.11     |
|                                                                                                                     |
| PetscLinearSolver                                                                                                   |
|   solve()                               3         0.9102      0.303389    0.9102      0.303389    38.69    38.69    |
|                                                                                                                     |
| ProjectVector                                                                                                       |
|   operator()                            2         0.0225      0.011261    0.1462      0.073105    0.96     6.22     |
|                                                                                                                     |
| System                                                                                                              |
|   assemble()                            3         0.1524      0.050816    0.4379      0.145950    6.48     18.61    |
|   project_vector()                      2         0.0032      0.001624    0.1548      0.077393    0.14     6.58     |
|                                                                                                                     |
| XdrIO                                                                                                               |
|   read()                                1         0.0003      0.000324    0.0004      0.000364    0.01     0.02     |
 ---------------------------------------------------------------------------------------------------------------------
| Totals:                                 144816    2.3523                                          100.00            |
 ---------------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example miscellaneous_ex5:
*  mpirun -np 4 example-devel  -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_ex5'

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

Hosted By:
SourceForge.net Logo