The source file subdomains_ex1.C with comments:

Subdomains Example 1 - Solving on a Subdomain



This example builds on the example 4 by showing what to do in order to solve an equation only on a subdomain.



C++ include files that we need
        #include <iostream>
        #include <algorithm>
        #include <math.h>
        
Basic include file needed for the mesh functionality.
        #include "libmesh/libmesh.h"
        #include "libmesh/mesh.h"
        #include "libmesh/mesh_generation.h"
        #include "libmesh/exodusII_io.h"
        #include "libmesh/gmv_io.h"
        #include "libmesh/gnuplot_io.h"
        #include "libmesh/linear_implicit_system.h"
        #include "libmesh/equation_systems.h"
        
Define the Finite Element object.
        #include "libmesh/fe.h"
        
Define Gauss quadrature rules.
        #include "libmesh/quadrature_gauss.h"
        
Define the DofMap, which handles degree of freedom indexing.
        #include "libmesh/dof_map.h"
        
Define useful datatypes for finite element matrix and vector components.
        #include "libmesh/sparse_matrix.h"
        #include "libmesh/numeric_vector.h"
        #include "libmesh/dense_matrix.h"
        #include "libmesh/dense_vector.h"
        
Define the PerfLog, a performance logging utility. It is useful for timing events in a code and giving you an idea where bottlenecks lie.
        #include "libmesh/perf_log.h"
        
The definition of a geometric element
        #include "libmesh/elem.h"
        
        #include "libmesh/mesh_refinement.h"
        
Classes needed for subdomain computation.
        #include "libmesh/system_subset_by_subdomain.h"
        
        #include "libmesh/string_to_enum.h"
        #include "libmesh/getpot.h"
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
        
        
Function prototype. This is the function that will assemble the linear system for our Poisson problem. Note that the function will take the \p EquationSystems object and the name of the system we are assembling as input. From the \p EquationSystems object we have acess to the \p Mesh and other objects we might need.
        void assemble_poisson(EquationSystems& es,
                              const std::string& system_name);
        
Exact solution function prototype.
        Real exact_solution (const Real x,
                             const Real y = 0.,
                             const Real z = 0.);
        
Begin the main program.
        int main (int argc, char** argv)
        {
Initialize libMesh and any dependent libaries, like in example 2.
          LibMeshInit init (argc, argv);
        
Only our PETSc interface currently supports solves restricted to subdomains
          libmesh_example_assert(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
        
Skip adaptive examples on a non-adaptive libMesh build
        #ifndef LIBMESH_ENABLE_AMR
          libmesh_example_assert(false, "--enable-amr");
        #else
        
Declare a performance log for the main program PerfLog perf_main("Main Program");

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] << " -d 2(3)" << " -n 15"
                          << std::endl;
        
This handy function will print the file name, line number, and then abort. Currrently the library does not use C++ exception handling.
              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 problem dimension from command line. Use int instead of unsigned since the GetPot overload is ambiguous otherwise.
          int dim = 2;
          if ( command_line.search(1, "-d") )
            dim = command_line.next(dim);
        
Skip higher-dimensional examples on a lower-dimensional libMesh build
          libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");
        
Create a mesh with user-defined dimension. Read number of elements from command line
          int ps = 15;
          if ( command_line.search(1, "-n") )
            ps = command_line.next(ps);
        
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 discontinuous basis.
          if ((family == "MONOMIAL") || (family == "XYZ"))
            {
              if (libMesh::processor_id() == 0)
                std::cerr << "ex4 currently requires a C^0 (or higher) FE basis." << std::endl;
              libmesh_error();
            }
        
Create a mesh, with dimension to be overridden later, on the default MPI communicator.
          Mesh mesh(init.comm());
        
Use the MeshTools::Generation mesh generator to create a uniform grid on the square [-1,1]^D. We instruct the mesh generator to build a mesh of 8x8 \p Quad9 elements in 2D, or \p Hex27 elements in 3D. Building these higher-order elements allows us to use higher-order approximation, as in example 3.

          Real halfwidth = dim > 1 ? 1. : 0.;
          Real halfheight = dim > 2 ? 1. : 0.;
        
          if ((family == "LAGRANGE") && (order == "FIRST"))
            {
No reason to use high-order geometric elements if we are solving with low-order finite elements.
              MeshTools::Generation::build_cube (mesh,
                                                 ps,
        					 (dim>1) ? ps : 0,
        					 (dim>2) ? ps : 0,
                                                 -1., 1.,
                                                 -halfwidth, halfwidth,
                                                 -halfheight, halfheight,
                                                 (dim==1)    ? EDGE2 :
                                                 ((dim == 2) ? QUAD4 : HEX8));
            }
        
          else
            {
              MeshTools::Generation::build_cube (mesh,
        					 ps,
        					 (dim>1) ? ps : 0,
        					 (dim>2) ? ps : 0,
                                                 -1., 1.,
                                                 -halfwidth, halfwidth,
                                                 -halfheight, halfheight,
                                                 (dim==1)    ? EDGE3 :
                                                 ((dim == 2) ? QUAD9 : HEX27));
            }
        
        
To demonstate solving on a subdomain, we will solve only on the interior of a circle (ball in 3d) with radius 0.8. So show that this also works well on locally refined meshes, we refine once all elements that are located on the boundary of this circle (or ball).
          {
A MeshRefinement object is needed to refine meshes.
            MeshRefinement meshRefinement(mesh);
        
Loop over all elements.
            MeshBase::element_iterator       elem_it  = mesh.elements_begin();
            const MeshBase::element_iterator elem_end = mesh.elements_end();
            for (; elem_it != elem_end; ++elem_it)
              {
        	Elem* elem = *elem_it;
        	if(elem->active())
        	  {
Just check whether the current element has at least one node inside and one node outside the circle.
                    bool node_in = false;
        	    bool node_out = false;
        	    for(unsigned int i=0; i<elem->n_nodes(); i++)
        	      {
        		double d = elem->point(i).size();
        		if(d<0.8)
        		  {
        		    node_in = true;
        		  }
        		else
        		  {
        		    node_out = true;
        		  }
        	      }
        	    if(node_in && node_out)
        	      {
        		elem->set_refinement_flag(Elem::REFINE);
        	      }
        	    else
        	      {
        		elem->set_refinement_flag(Elem::DO_NOTHING);
        	      }
        	  }
        	else
        	  {
        	    elem->set_refinement_flag(Elem::INACTIVE);
        	  }
              }
        
Now actually refine.
            meshRefinement.refine_elements();
          }
        
Print information about the mesh to the screen.
          mesh.print_info();
        
Now set the subdomain_id of all elements whose centroid is inside the circle to 1.
          {
Loop over all elements.
            MeshBase::element_iterator       elem_it  = mesh.elements_begin();
            const MeshBase::element_iterator elem_end = mesh.elements_end();
            for (; elem_it != elem_end; ++elem_it)
              {
        	Elem* elem = *elem_it;
        	double d = elem->centroid().size();
        	if(d<0.8)
        	  {
        	    elem->subdomain_id() = 1;
        	  }
              }
          }
        
Create an equation systems object.
          EquationSystems equation_systems (mesh);
        
Declare the system and its variables. Create a system named "Poisson"
          LinearImplicitSystem& system =
            equation_systems.add_system<LinearImplicitSystem> ("Poisson");
        
        
Add the variable "u" to "Poisson". "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 matrix assembly function.
          system.attach_assemble_function (assemble_poisson);
        
Initialize the data structures for the equation system.
          equation_systems.init();
        
Print information about the system to the screen.
          equation_systems.print_info();
          mesh.print_info();
        
Restrict solves to those elements that have subdomain_id set to 1.
          std::set<subdomain_id_type> id_list;
          id_list.insert(1);
          SystemSubsetBySubdomain::SubdomainSelectionByList selection(id_list);
          SystemSubsetBySubdomain subset(system,selection);
          system.restrict_solve_to(&subset,SUBSET_ZERO);
        
Note that using \p SUBSET_ZERO will cause all dofs outside the subdomain to be cleared. This will, however, cause some hanging nodes outside the subdomain to have inconsistent values.

Solve the system "Poisson", just like example 2.
          equation_systems.get_system("Poisson").solve();
        
After solving the system write the solution to a GMV-formatted plot file.
          if(dim == 1)
          {
            GnuPlotIO plot(mesh,"Subdomains Example 1, 1D",GnuPlotIO::GRID_ON);
            plot.write_equation_systems("gnuplot_script",equation_systems);
          }
          else
          {
            GMVIO (mesh).write_equation_systems ((dim == 3) ?
              "out_3.gmv" : "out_2.gmv",equation_systems);
        #ifdef LIBMESH_HAVE_EXODUS_API
            ExodusII_IO (mesh).write_equation_systems ((dim == 3) ?
              "out_3.e" : "out_2.e",equation_systems);
        #endif // #ifdef LIBMESH_HAVE_EXODUS_API
          }
        
        #endif // #ifndef LIBMESH_ENABLE_AMR
        
All done.
          return 0;
        }
        
        
        
        






We now define the matrix assembly function for the Poisson system. We need to first compute element matrices and right-hand sides, and then take into account the boundary conditions, which will be handled via a penalty method.
        void assemble_poisson(EquationSystems& es,
                              const std::string& system_name)
        {
It is a good idea to make sure we are assembling the proper system.
          libmesh_assert_equal_to (system_name, "Poisson");
        
        
Declare a performance log. Give it a descriptive string to identify what part of the code we are logging, since there may be many PerfLogs in an application.
          PerfLog perf_log ("Matrix Assembly");
        
Get a constant reference to the mesh object.
          const MeshBase& mesh = es.get_mesh();
        
The dimension that we are running
          const unsigned int dim = mesh.mesh_dimension();
        
Get a reference to the LinearImplicitSystem we are solving
          LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Poisson");
        
A reference to the \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 physical XY locations of the quadrature points on the element. These might be useful for evaluating spatially varying material properties at the quadrature points.
          const std::vector<Point>& q_point = fe->get_xyz();
        
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 element matrix and right-hand-side vector contribution. Following basic finite element terminology we will denote these "Ke" and "Fe". More detail is in example 3.
          DenseMatrix<Number> Ke;
          DenseVector<Number> Fe;
        
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 the element matrix and right-hand-side contribution. See example 3 for a discussion of the element iterators.
          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;
        
Elements with subdomain_id other than 1 are not in the active subdomain. We don't assemble anything for them.
              if(elem->subdomain_id()==1)
        	{
Start logging the shape function initialization. This is done through a simple function call with the name of the event to log.
                  perf_log.push("elem init");
        
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 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. 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());
        
        	  Fe.resize (dof_indices.size());
        
Stop logging the shape function initialization. If you forget to stop logging an event the PerfLog object will probably catch the error and abort.
                  perf_log.pop("elem init");
        
Now we will build the element matrix. This involves a double loop to integrate the test funcions (i) against the trial functions (j).

We have split the numeric integration into two loops so that we can log the matrix and right-hand-side computation seperately.

Now start logging the element matrix computation
                  perf_log.push ("Ke");
        
        	  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        	    for (unsigned int i=0; i<phi.size(); i++)
        	      for (unsigned int j=0; j<phi.size(); j++)
        		Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
        
        
Stop logging the matrix computation
                  perf_log.pop ("Ke");
        
Now we build the element right-hand-side contribution. This involves a single loop in which we integrate the "forcing function" in the PDE against the test functions.

Start logging the right-hand-side computation
                  perf_log.push ("Fe");
        
        	  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
        	    {
fxy is the forcing function for the Poisson equation. In this case we set fxy to be a finite difference Laplacian approximation to the (known) exact solution.

We will use the second-order accurate FD Laplacian approximation, which in 2D on a structured grid is

u_xx + u_yy = (u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) + -4*u(i,j))/h^2

Since the value of the forcing function depends only on the location of the quadrature point (q_point[qp]) we will compute it here, outside of the i-loop
                      const Real x = q_point[qp](0);
        #if LIBMESH_DIM > 1
        	      const Real y = q_point[qp](1);
        #else
        	      const Real y = 0.;
        #endif
        #if LIBMESH_DIM > 2
        	      const Real z = q_point[qp](2);
        #else
        	      const Real z = 0.;
        #endif
        	      const Real eps = 1.e-3;
        
        	      const Real uxx = (exact_solution(x-eps,y,z) +
        				exact_solution(x+eps,y,z) +
        				-2.*exact_solution(x,y,z))/eps/eps;
        
        	      const Real uyy = (exact_solution(x,y-eps,z) +
        				exact_solution(x,y+eps,z) +
        				-2.*exact_solution(x,y,z))/eps/eps;
        
        	      const Real uzz = (exact_solution(x,y,z-eps) +
        				exact_solution(x,y,z+eps) +
        				-2.*exact_solution(x,y,z))/eps/eps;
        
        	      Real fxy;
        	      if(dim==1)
        		{
In 1D, compute the rhs by differentiating the exact solution twice.
                          const Real pi = libMesh::pi;
        		  fxy = (0.25*pi*pi)*sin(.5*pi*x);
        		}
        	      else
        		{
        		  fxy = - (uxx + uyy + ((dim==2) ? 0. : uzz));
        		}
        
Add the RHS contribution
                      for (unsigned int i=0; i<phi.size(); i++)
        		Fe(i) += JxW[qp]*fxy*phi[i][qp];
        	    }
        
Stop logging the right-hand-side computation
                  perf_log.pop ("Fe");
        
At this point the interior element integration has been completed. However, we have not yet addressed boundary conditions. For this example we will only consider simple Dirichlet boundary conditions imposed via the penalty method. This is discussed at length in example 3.
                  {
        
Start logging the boundary condition computation
                    perf_log.push ("BCs");
        
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. If there is a neighbor, check that neighbor's subdomain_id; if that is different from 1, the side is also located on the boundary.
                    for (unsigned int side=0; side<elem->n_sides(); side++)
        	      if ((elem->neighbor(side) == NULL) ||
        		  (elem->neighbor(side)->subdomain_id()!=1))
        		{
        
The penalty value. \frac{1}{\epsilon} in the discussion above.
                          const Real penalty = 1.e10;
        
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();
        
The XYZ locations (in physical space) of the quadrature points on the face. This is where we will interpolate the boundary value function.
                          const std::vector<Point >& qface_point = fe_face->get_xyz();
        
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++)
        		    {
The location on the boundary of the current face quadrature point.
                              const Real xf = qface_point[qp](0);
        #if LIBMESH_DIM > 1
        		      const Real yf = qface_point[qp](1);
        #else
        		      const Real yf = 0.;
        #endif
        #if LIBMESH_DIM > 2
        		      const Real zf = qface_point[qp](2);
        #else
        		      const Real zf = 0.;
        #endif
        
        
The boundary value.
                              const Real value = exact_solution(xf, yf, zf);
        
Matrix contribution of the L2 projection.
                              for (unsigned int i=0; i<phi_face.size(); i++)
        			for (unsigned int j=0; j<phi_face.size(); j++)
        			  Ke(i,j) += JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp];
        
Right-hand-side contribution of the L2 projection.
                              for (unsigned int i=0; i<phi_face.size(); i++)
        			Fe(i) += JxW_face[qp]*penalty*value*phi_face[i][qp];
        		    }
        		}
        
        
Stop logging the boundary condition computation
                    perf_log.pop ("BCs");
        	  }
        
If this assembly program were to be used on an adaptive mesh, we would have to apply any hanging node constraint equations
                  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
        
The element 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. Start logging the insertion of the local (element) matrix and vector into the global matrix and vector
                  perf_log.push ("matrix insertion");
        
        	  system.matrix->add_matrix (Ke, dof_indices);
        	  system.rhs->add_vector    (Fe, dof_indices);
        
Start logging the insertion of the local (element) matrix and vector into the global matrix and vector
                  perf_log.pop ("matrix insertion");
        	}
            }
        
That's it. We don't need to do anything else to the PerfLog. When it goes out of scope (at this function return) it will print its log to the screen. Pretty easy, huh?
        }



The source file subdomains_ex1.C without comments:

 
  
  
  #include <iostream>
  #include <algorithm>
  #include <math.h>
  
  #include "libmesh/libmesh.h"
  #include "libmesh/mesh.h"
  #include "libmesh/mesh_generation.h"
  #include "libmesh/exodusII_io.h"
  #include "libmesh/gmv_io.h"
  #include "libmesh/gnuplot_io.h"
  #include "libmesh/linear_implicit_system.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/perf_log.h"
  
  #include "libmesh/elem.h"
  
  #include "libmesh/mesh_refinement.h"
  
  #include "libmesh/system_subset_by_subdomain.h"
  
  #include "libmesh/string_to_enum.h"
  #include "libmesh/getpot.h"
  
  using namespace libMesh;
  
  
  
  void assemble_poisson(EquationSystems& es,
                        const std::string& system_name);
  
  Real exact_solution (const Real x,
                       const Real y = 0.,
                       const Real z = 0.);
  
  int main (int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
  
    libmesh_example_assert(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
  
  #ifndef LIBMESH_ENABLE_AMR
    libmesh_example_assert(false, "--enable-amr");
  #else
  
  
    GetPot command_line (argc, argv);
  
    if (argc < 3)
      {
        if (libMesh::processor_id() == 0)
          std::cerr << "Usage:\n"
                    <<"\t " << argv[0] << " -d 2(3)" << " -n 15"
                    << 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 dim = 2;
    if ( command_line.search(1, "-d") )
      dim = command_line.next(dim);
  
    libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");
  
    int ps = 15;
    if ( command_line.search(1, "-n") )
      ps = command_line.next(ps);
  
    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"))
      {
        if (libMesh::processor_id() == 0)
          std::cerr << "ex4 currently requires a C^0 (or higher) FE basis." << std::endl;
        libmesh_error();
      }
  
    Mesh mesh(init.comm());
  
  
    Real halfwidth = dim > 1 ? 1. : 0.;
    Real halfheight = dim > 2 ? 1. : 0.;
  
    if ((family == "LAGRANGE") && (order == "FIRST"))
      {
        MeshTools::Generation::build_cube (mesh,
                                           ps,
  					 (dim>1) ? ps : 0,
  					 (dim>2) ? ps : 0,
                                           -1., 1.,
                                           -halfwidth, halfwidth,
                                           -halfheight, halfheight,
                                           (dim==1)    ? EDGE2 :
                                           ((dim == 2) ? QUAD4 : HEX8));
      }
  
    else
      {
        MeshTools::Generation::build_cube (mesh,
  					 ps,
  					 (dim>1) ? ps : 0,
  					 (dim>2) ? ps : 0,
                                           -1., 1.,
                                           -halfwidth, halfwidth,
                                           -halfheight, halfheight,
                                           (dim==1)    ? EDGE3 :
                                           ((dim == 2) ? QUAD9 : HEX27));
      }
  
  
    {
      MeshRefinement meshRefinement(mesh);
  
      MeshBase::element_iterator       elem_it  = mesh.elements_begin();
      const MeshBase::element_iterator elem_end = mesh.elements_end();
      for (; elem_it != elem_end; ++elem_it)
        {
  	Elem* elem = *elem_it;
  	if(elem->active())
  	  {
  	    bool node_in = false;
  	    bool node_out = false;
  	    for(unsigned int i=0; i<elem->n_nodes(); i++)
  	      {
  		double d = elem->point(i).size();
  		if(d<0.8)
  		  {
  		    node_in = true;
  		  }
  		else
  		  {
  		    node_out = true;
  		  }
  	      }
  	    if(node_in && node_out)
  	      {
  		elem->set_refinement_flag(Elem::REFINE);
  	      }
  	    else
  	      {
  		elem->set_refinement_flag(Elem::DO_NOTHING);
  	      }
  	  }
  	else
  	  {
  	    elem->set_refinement_flag(Elem::INACTIVE);
  	  }
        }
  
      meshRefinement.refine_elements();
    }
  
    mesh.print_info();
  
    {
      MeshBase::element_iterator       elem_it  = mesh.elements_begin();
      const MeshBase::element_iterator elem_end = mesh.elements_end();
      for (; elem_it != elem_end; ++elem_it)
        {
  	Elem* elem = *elem_it;
  	double d = elem->centroid().size();
  	if(d<0.8)
  	  {
  	    elem->subdomain_id() = 1;
  	  }
        }
    }
  
    EquationSystems equation_systems (mesh);
  
    LinearImplicitSystem& system =
      equation_systems.add_system<LinearImplicitSystem> ("Poisson");
  
  
    system.add_variable("u",
                        Utility::string_to_enum<Order>   (order),
                        Utility::string_to_enum<FEFamily>(family));
  
    system.attach_assemble_function (assemble_poisson);
  
    equation_systems.init();
  
    equation_systems.print_info();
    mesh.print_info();
  
    std::set<subdomain_id_type> id_list;
    id_list.insert(1);
    SystemSubsetBySubdomain::SubdomainSelectionByList selection(id_list);
    SystemSubsetBySubdomain subset(system,selection);
    system.restrict_solve_to(&subset,SUBSET_ZERO);
  
  
    equation_systems.get_system("Poisson").solve();
  
    if(dim == 1)
    {
      GnuPlotIO plot(mesh,"Subdomains Example 1, 1D",GnuPlotIO::GRID_ON);
      plot.write_equation_systems("gnuplot_script",equation_systems);
    }
    else
    {
      GMVIO (mesh).write_equation_systems ((dim == 3) ?
        "out_3.gmv" : "out_2.gmv",equation_systems);
  #ifdef LIBMESH_HAVE_EXODUS_API
      ExodusII_IO (mesh).write_equation_systems ((dim == 3) ?
        "out_3.e" : "out_2.e",equation_systems);
  #endif // #ifdef LIBMESH_HAVE_EXODUS_API
    }
  
  #endif // #ifndef LIBMESH_ENABLE_AMR
  
    return 0;
  }
  
  
  
  
  void assemble_poisson(EquationSystems& es,
                        const std::string& system_name)
  {
    libmesh_assert_equal_to (system_name, "Poisson");
  
  
    PerfLog perf_log ("Matrix Assembly");
  
    const MeshBase& mesh = es.get_mesh();
  
    const unsigned int dim = mesh.mesh_dimension();
  
    LinearImplicitSystem& system = es.get_system<LinearImplicitSystem>("Poisson");
  
    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<Point>& q_point = fe->get_xyz();
  
    const std::vector<std::vector<Real> >& phi = fe->get_phi();
  
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
    DenseMatrix<Number> Ke;
    DenseVector<Number> Fe;
  
    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;
  
        if(elem->subdomain_id()==1)
  	{
  	  perf_log.push("elem init");
  
  	  dof_map.dof_indices (elem, dof_indices);
  
  	  fe->reinit (elem);
  
  	  Ke.resize (dof_indices.size(),
  		     dof_indices.size());
  
  	  Fe.resize (dof_indices.size());
  
  	  perf_log.pop("elem init");
  
  	  perf_log.push ("Ke");
  
  	  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
  	    for (unsigned int i=0; i<phi.size(); i++)
  	      for (unsigned int j=0; j<phi.size(); j++)
  		Ke(i,j) += JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
  
  
  	  perf_log.pop ("Ke");
  
  	  perf_log.push ("Fe");
  
  	  for (unsigned int qp=0; qp<qrule.n_points(); qp++)
  	    {
  	      const Real x = q_point[qp](0);
  #if LIBMESH_DIM > 1
  	      const Real y = q_point[qp](1);
  #else
  	      const Real y = 0.;
  #endif
  #if LIBMESH_DIM > 2
  	      const Real z = q_point[qp](2);
  #else
  	      const Real z = 0.;
  #endif
  	      const Real eps = 1.e-3;
  
  	      const Real uxx = (exact_solution(x-eps,y,z) +
  				exact_solution(x+eps,y,z) +
  				-2.*exact_solution(x,y,z))/eps/eps;
  
  	      const Real uyy = (exact_solution(x,y-eps,z) +
  				exact_solution(x,y+eps,z) +
  				-2.*exact_solution(x,y,z))/eps/eps;
  
  	      const Real uzz = (exact_solution(x,y,z-eps) +
  				exact_solution(x,y,z+eps) +
  				-2.*exact_solution(x,y,z))/eps/eps;
  
  	      Real fxy;
  	      if(dim==1)
  		{
  		  const Real pi = libMesh::pi;
  		  fxy = (0.25*pi*pi)*sin(.5*pi*x);
  		}
  	      else
  		{
  		  fxy = - (uxx + uyy + ((dim==2) ? 0. : uzz));
  		}
  
  	      for (unsigned int i=0; i<phi.size(); i++)
  		Fe(i) += JxW[qp]*fxy*phi[i][qp];
  	    }
  
  	  perf_log.pop ("Fe");
  
  	  {
  
  	    perf_log.push ("BCs");
  
  	    for (unsigned int side=0; side<elem->n_sides(); side++)
  	      if ((elem->neighbor(side) == NULL) ||
  		  (elem->neighbor(side)->subdomain_id()!=1))
  		{
  
  		  const Real penalty = 1.e10;
  
  		  const std::vector<std::vector<Real> >&  phi_face = fe_face->get_phi();
  
  		  const std::vector<Real>& JxW_face = fe_face->get_JxW();
  
  		  const std::vector<Point >& qface_point = fe_face->get_xyz();
  
  		  fe_face->reinit(elem, side);
  
  		  for (unsigned int qp=0; qp<qface.n_points(); qp++)
  		    {
  		      const Real xf = qface_point[qp](0);
  #if LIBMESH_DIM > 1
  		      const Real yf = qface_point[qp](1);
  #else
  		      const Real yf = 0.;
  #endif
  #if LIBMESH_DIM > 2
  		      const Real zf = qface_point[qp](2);
  #else
  		      const Real zf = 0.;
  #endif
  
  
  		      const Real value = exact_solution(xf, yf, zf);
  
  		      for (unsigned int i=0; i<phi_face.size(); i++)
  			for (unsigned int j=0; j<phi_face.size(); j++)
  			  Ke(i,j) += JxW_face[qp]*penalty*phi_face[i][qp]*phi_face[j][qp];
  
  		      for (unsigned int i=0; i<phi_face.size(); i++)
  			Fe(i) += JxW_face[qp]*penalty*value*phi_face[i][qp];
  		    }
  		}
  
  
  	    perf_log.pop ("BCs");
  	  }
  
  	  dof_map.constrain_element_matrix_and_vector (Ke, Fe, dof_indices);
  
  	  perf_log.push ("matrix insertion");
  
  	  system.matrix->add_matrix (Ke, dof_indices);
  	  system.rhs->add_vector    (Fe, dof_indices);
  
  	  perf_log.pop ("matrix insertion");
  	}
      }
  
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/subdomains/subdomains_ex1'
***************************************************************
* Running Example subdomains_ex1:
*  mpirun -np 4 example-devel -d 1 -n 40 -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/subdomains/subdomains_ex1/.libs/lt-example-devel -d 1 -n 40 -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc

 Mesh Information:
  mesh_dimension()=1
  spatial_dimension()=3
  n_nodes()=43
    n_local_nodes()=12
  n_elem()=44
    n_local_elem()=12
    n_active_elem()=42
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "Poisson"
    Type "LinearImplicit"
    Variables="u" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="FIRST", "THIRD" 
    n_dofs()=43
    n_local_dofs()=12
    n_constrained_dofs()=0
    n_local_constrained_dofs()=0
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 2.88372
      Average Off-Processor Bandwidth <= 0.139535
      Maximum  On-Processor Bandwidth <= 3
      Maximum Off-Processor Bandwidth <= 1
    DofMap Constraints
      Number of DoF Constraints = 0
      Number of Node Constraints = 0

 Mesh Information:
  mesh_dimension()=1
  spatial_dimension()=3
  n_nodes()=43
    n_local_nodes()=12
  n_elem()=44
    n_local_elem()=12
    n_active_elem()=42
  n_subdomains()=2
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0


 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:52:31 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'                                               |
 -------------------------------------------------------------------------------------------------------------------
 -----------------------------------------------------------------------------------------------------------
| Matrix Assembly Performance: Alive time=0.001667, Active time=0.001447                                    |
 -----------------------------------------------------------------------------------------------------------
| 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   |
|-----------------------------------------------------------------------------------------------------------|
|                                                                                                           |
| BCs                           7         0.0001      0.000009    0.0001      0.000009    4.15     4.15     |
| Fe                            7         0.0000      0.000006    0.0000      0.000006    3.04     3.04     |
| Ke                            7         0.0000      0.000001    0.0000      0.000001    0.48     0.48     |
| elem init                     7         0.0013      0.000186    0.0013      0.000186    89.84    89.84    |
| matrix insertion              7         0.0000      0.000005    0.0000      0.000005    2.49     2.49     |
 -----------------------------------------------------------------------------------------------------------
| Totals:                       35        0.0014                                          100.00            |
 -----------------------------------------------------------------------------------------------------------

 ----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=0.212075, Active time=0.143039                                                 |
 ----------------------------------------------------------------------------------------------------------------
| 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.0001      0.000074    0.0001      0.000085    0.05     0.06     |
|   build_sparsity()                 1         0.0001      0.000127    0.0004      0.000397    0.09     0.28     |
|   create_dof_constraints()         1         0.0000      0.000002    0.0000      0.000002    0.00     0.00     |
|   distribute_dofs()                1         0.0002      0.000245    0.0097      0.009703    0.17     6.78     |
|   dof_indices()                    38        0.0001      0.000003    0.0001      0.000003    0.09     0.09     |
|   prepare_send_list()              1         0.0000      0.000002    0.0000      0.000002    0.00     0.00     |
|   reinit()                         1         0.0009      0.000895    0.0009      0.000895    0.63     0.63     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          1         0.0001      0.000132    0.0041      0.004078    0.09     2.85     |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        8         0.0000      0.000003    0.0000      0.000003    0.02     0.02     |
|   init_shape_functions()           2         0.0000      0.000008    0.0000      0.000008    0.01     0.01     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             8         0.0000      0.000003    0.0000      0.000003    0.02     0.02     |
|   compute_face_map()               1         0.0000      0.000006    0.0000      0.000006    0.00     0.00     |
|   init_face_shape_functions()      1         0.0000      0.000006    0.0000      0.000006    0.00     0.00     |
|   init_reference_to_physical_map() 2         0.0000      0.000008    0.0000      0.000008    0.01     0.01     |
|                                                                                                                |
| GnuPlotIO                                                                                                      |
|   write_nodal_data()               1         0.0006      0.000625    0.0006      0.000625    0.44     0.44     |
|                                                                                                                |
| LocationMap                                                                                                    |
|   find()                           4         0.0000      0.000002    0.0000      0.000002    0.01     0.01     |
|   init()                           1         0.0004      0.000382    0.0004      0.000382    0.27     0.27     |
|                                                                                                                |
| Mesh                                                                                                           |
|   find_neighbors()                 2         0.0003      0.000172    0.0051      0.002543    0.24     3.56     |
|   renumber_nodes_and_elem()        4         0.0000      0.000009    0.0000      0.000009    0.03     0.03     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   compute_hilbert_indices()        3         0.0006      0.000203    0.0006      0.000203    0.43     0.43     |
|   find_global_indices()            3         0.0002      0.000069    0.0069      0.002286    0.14     4.79     |
|   parallel_sort()                  3         0.0016      0.000539    0.0036      0.001203    1.13     2.52     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         1         0.0001      0.000069    0.0062      0.006247    0.05     4.37     |
|                                                                                                                |
| MeshRefinement                                                                                                 |
|   _refine_elements()               1         0.0001      0.000057    0.0002      0.000161    0.04     0.11     |
|   add_point()                      4         0.0000      0.000004    0.0000      0.000007    0.01     0.02     |
|   make_flags_parallel_consistent() 1         0.0001      0.000098    0.0080      0.008026    0.07     5.61     |
|   make_refinement_compatible()     2         0.0000      0.000014    0.0006      0.000287    0.02     0.40     |
|                                                                                                                |
| MeshTools::Generation                                                                                          |
|   build_cube()                     1         0.0001      0.000098    0.0001      0.000098    0.07     0.07     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      2         0.0066      0.003286    0.0184      0.009180    4.59     12.84    |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      13        0.0024      0.000185    0.0032      0.000248    1.68     2.25     |
|   max(bool)                        5         0.0009      0.000173    0.0009      0.000173    0.60     0.60     |
|   max(scalar)                      182       0.0255      0.000140    0.0255      0.000140    17.86    17.86    |
|   max(vector)                      43        0.0066      0.000154    0.0254      0.000591    4.62     17.76    |
|   min(bool)                        218       0.0327      0.000150    0.0327      0.000150    22.87    22.87    |
|   min(scalar)                      175       0.0279      0.000159    0.0279      0.000159    19.51    19.51    |
|   min(vector)                      43        0.0066      0.000153    0.0260      0.000604    4.59     18.16    |
|   probe()                          63        0.0047      0.000074    0.0047      0.000074    3.28     3.28     |
|   receive()                        63        0.0002      0.000003    0.0049      0.000078    0.14     3.42     |
|   send()                           63        0.0001      0.000002    0.0001      0.000002    0.08     0.08     |
|   send_receive()                   66        0.0003      0.000004    0.0053      0.000080    0.20     3.70     |
|   sum()                            23        0.0016      0.000068    0.0033      0.000142    1.09     2.28     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           63        0.0001      0.000001    0.0001      0.000001    0.06     0.06     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         2         0.0002      0.000097    0.0089      0.004469    0.14     6.25     |
|   set_parent_processor_ids()       2         0.0000      0.000022    0.0000      0.000022    0.03     0.03     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          1         0.0190      0.019050    0.0190      0.019050    13.32    13.32    |
|                                                                                                                |
| System                                                                                                         |
|   assemble()                       1         0.0017      0.001737    0.0019      0.001856    1.21     1.30     |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            1126      0.1430                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example subdomains_ex1:
*  mpirun -np 4 example-devel -d 1 -n 40 -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
***************************************************************
* Running Example subdomains_ex1:
*  mpirun -np 4 example-devel -d 2 -n 30 -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/subdomains/subdomains_ex1/.libs/lt-example-devel -d 2 -n 30 -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()=1329
    n_local_nodes()=349
  n_elem()=1268
    n_local_elem()=315
    n_active_elem()=1176
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "Poisson"
    Type "LinearImplicit"
    Variables="u" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="FIRST", "THIRD" 
    n_dofs()=1329
    n_local_dofs()=349
    n_constrained_dofs()=184
    n_local_constrained_dofs()=42
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 9.03988
      Average Off-Processor Bandwidth <= 0.337096
      Maximum  On-Processor Bandwidth <= 15
      Maximum Off-Processor Bandwidth <= 6
    DofMap Constraints
      Number of DoF Constraints = 184
      Average DoF Constraint Length= 2
      Number of Node Constraints = 368
      Maximum Node Constraint Length= 3
      Average Node Constraint Length= 2.5

 Mesh Information:
  mesh_dimension()=2
  spatial_dimension()=3
  n_nodes()=1329
    n_local_nodes()=349
  n_elem()=1268
    n_local_elem()=315
    n_active_elem()=1176
  n_subdomains()=2
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0


 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:52:31 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'                                               |
 -------------------------------------------------------------------------------------------------------------------
 -----------------------------------------------------------------------------------------------------------
| Matrix Assembly Performance: Alive time=0.015978, Active time=0.007801                                    |
 -----------------------------------------------------------------------------------------------------------
| 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   |
|-----------------------------------------------------------------------------------------------------------|
|                                                                                                           |
| BCs                           163       0.0021      0.000013    0.0021      0.000013    26.34    26.34    |
| Fe                            163       0.0013      0.000008    0.0013      0.000008    16.15    16.15    |
| Ke                            163       0.0008      0.000005    0.0008      0.000005    10.40    10.40    |
| elem init                     163       0.0032      0.000020    0.0032      0.000020    40.92    40.92    |
| matrix insertion              163       0.0005      0.000003    0.0005      0.000003    6.19     6.19     |
 -----------------------------------------------------------------------------------------------------------
| Totals:                       815       0.0078                                          100.00            |
 -----------------------------------------------------------------------------------------------------------

 ----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=0.214073, Active time=0.170199                                                 |
 ----------------------------------------------------------------------------------------------------------------
| 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.0016      0.001580    0.0020      0.001981    0.93     1.16     |
|   build_constraint_matrix()        163       0.0004      0.000002    0.0004      0.000002    0.22     0.22     |
|   build_sparsity()                 1         0.0015      0.001474    0.0034      0.003370    0.87     1.98     |
|   cnstrn_elem_mat_vec()            163       0.0065      0.000040    0.0065      0.000040    3.81     3.81     |
|   create_dof_constraints()         1         0.0039      0.003922    0.0077      0.007718    2.30     4.53     |
|   distribute_dofs()                1         0.0033      0.003334    0.0091      0.009094    1.96     5.34     |
|   dof_indices()                    2021      0.0088      0.000004    0.0088      0.000004    5.16     5.16     |
|   prepare_send_list()              1         0.0000      0.000007    0.0000      0.000007    0.00     0.00     |
|   reinit()                         1         0.0054      0.005430    0.0054      0.005430    3.19     3.19     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          2         0.0017      0.000842    0.0053      0.002649    0.99     3.11     |
|                                                                                                                |
| ExodusII_IO                                                                                                    |
|   write_nodal_data()               1         0.0420      0.042039    0.0420      0.042039    24.70    24.70    |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        209       0.0010      0.000005    0.0010      0.000005    0.60     0.60     |
|   init_shape_functions()           47        0.0001      0.000002    0.0001      0.000002    0.06     0.06     |
|   inverse_map()                    1058      0.0019      0.000002    0.0019      0.000002    1.09     1.09     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             209       0.0007      0.000003    0.0007      0.000003    0.40     0.40     |
|   compute_face_map()               46        0.0003      0.000008    0.0009      0.000019    0.20     0.51     |
|   init_face_shape_functions()      1         0.0000      0.000012    0.0000      0.000012    0.01     0.01     |
|   init_reference_to_physical_map() 47        0.0003      0.000006    0.0003      0.000006    0.15     0.15     |
|                                                                                                                |
| GMVIO                                                                                                          |
|   write_nodal_data()               1         0.0103      0.010344    0.0104      0.010419    6.08     6.12     |
|                                                                                                                |
| LocationMap                                                                                                    |
|   find()                           1104      0.0014      0.000001    0.0014      0.000001    0.81     0.81     |
|   init()                           1         0.0004      0.000449    0.0004      0.000449    0.26     0.26     |
|                                                                                                                |
| Mesh                                                                                                           |
|   find_neighbors()                 2         0.0065      0.003235    0.0095      0.004745    3.80     5.58     |
|   renumber_nodes_and_elem()        4         0.0005      0.000127    0.0005      0.000127    0.30     0.30     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   compute_hilbert_indices()        3         0.0127      0.004227    0.0127      0.004227    7.45     7.45     |
|   find_global_indices()            3         0.0018      0.000608    0.0157      0.005234    1.07     9.23     |
|   parallel_sort()                  3         0.0006      0.000186    0.0008      0.000255    0.33     0.45     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         2         0.0001      0.000050    0.0580      0.029002    0.06     34.08    |
|                                                                                                                |
| MeshRefinement                                                                                                 |
|   _refine_elements()               1         0.0020      0.001962    0.0059      0.005896    1.15     3.46     |
|   add_point()                      1104      0.0020      0.000002    0.0035      0.000003    1.19     2.07     |
|   make_flags_parallel_consistent() 1         0.0006      0.000627    0.0009      0.000929    0.37     0.55     |
|   make_refinement_compatible()     2         0.0003      0.000161    0.0003      0.000166    0.19     0.20     |
|                                                                                                                |
| MeshTools::Generation                                                                                          |
|   build_cube()                     1         0.0007      0.000716    0.0007      0.000716    0.42     0.42     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      2         0.0220      0.011015    0.0329      0.016454    12.94    19.34    |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      13        0.0002      0.000012    0.0002      0.000017    0.09     0.13     |
|   max(bool)                        5         0.0003      0.000063    0.0003      0.000063    0.19     0.19     |
|   max(scalar)                      201       0.0021      0.000010    0.0021      0.000010    1.21     1.21     |
|   max(vector)                      47        0.0006      0.000014    0.0018      0.000038    0.37     1.05     |
|   min(bool)                        241       0.0020      0.000008    0.0020      0.000008    1.20     1.20     |
|   min(scalar)                      194       0.0063      0.000033    0.0063      0.000033    3.73     3.73     |
|   min(vector)                      47        0.0007      0.000014    0.0020      0.000043    0.39     1.19     |
|   probe()                          63        0.0002      0.000003    0.0002      0.000003    0.12     0.12     |
|   receive()                        63        0.0002      0.000003    0.0004      0.000007    0.13     0.25     |
|   send()                           63        0.0001      0.000002    0.0001      0.000002    0.08     0.08     |
|   send_receive()                   66        0.0003      0.000005    0.0009      0.000014    0.19     0.55     |
|   sum()                            26        0.0004      0.000015    0.0006      0.000023    0.23     0.35     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           63        0.0001      0.000001    0.0001      0.000001    0.05     0.05     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         2         0.0020      0.000981    0.0023      0.001128    1.15     1.33     |
|   set_parent_processor_ids()       2         0.0007      0.000361    0.0007      0.000361    0.42     0.42     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          1         0.0072      0.007205    0.0072      0.007205    4.23     4.23     |
|                                                                                                                |
| System                                                                                                         |
|   assemble()                       1         0.0054      0.005399    0.0162      0.016209    3.17     9.52     |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            7305      0.1702                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example subdomains_ex1:
*  mpirun -np 4 example-devel -d 2 -n 30 -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
***************************************************************
* Running Example subdomains_ex1:
*  mpirun -np 4 example-devel -d 3 -n 15 -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/subdomains/subdomains_ex1/.libs/lt-example-devel -d 3 -n 15 -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc

 Mesh Information:
  mesh_dimension()=3
  spatial_dimension()=3
  n_nodes()=10854
    n_local_nodes()=2973
  n_elem()=8767
    n_local_elem()=2192
    n_active_elem()=8093
  n_subdomains()=1
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0

 EquationSystems
  n_systems()=1
   System #0, "Poisson"
    Type "LinearImplicit"
    Variables="u" 
    Finite Element Types="LAGRANGE", "JACOBI_20_00" 
    Infinite Element Mapping="CARTESIAN" 
    Approximation Orders="FIRST", "THIRD" 
    n_dofs()=10854
    n_local_dofs()=2973
    n_constrained_dofs()=4068
    n_local_constrained_dofs()=1074
    n_vectors()=1
    n_matrices()=1
    DofMap Sparsity
      Average  On-Processor Bandwidth <= 28.3958
      Average Off-Processor Bandwidth <= 2.17652
      Maximum  On-Processor Bandwidth <= 78
      Maximum Off-Processor Bandwidth <= 49
    DofMap Constraints
      Number of DoF Constraints = 4068
      Average DoF Constraint Length= 2.66667
      Number of Node Constraints = 5428
      Maximum Node Constraint Length= 13
      Average Node Constraint Length= 6.24613

 Mesh Information:
  mesh_dimension()=3
  spatial_dimension()=3
  n_nodes()=10854
    n_local_nodes()=2973
  n_elem()=8767
    n_local_elem()=2192
    n_active_elem()=8093
  n_subdomains()=2
  n_partitions()=4
  n_processors()=4
  n_threads()=1
  processor_id()=0


 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:52:33 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'                                               |
 -------------------------------------------------------------------------------------------------------------------
 -----------------------------------------------------------------------------------------------------------
| Matrix Assembly Performance: Alive time=0.192523, Active time=0.162361                                    |
 -----------------------------------------------------------------------------------------------------------
| 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   |
|-----------------------------------------------------------------------------------------------------------|
|                                                                                                           |
| BCs                           745       0.0568      0.000076    0.0568      0.000076    35.01    35.01    |
| Fe                            745       0.0186      0.000025    0.0186      0.000025    11.45    11.45    |
| Ke                            745       0.0353      0.000047    0.0353      0.000047    21.75    21.75    |
| elem init                     745       0.0442      0.000059    0.0442      0.000059    27.24    27.24    |
| matrix insertion              745       0.0074      0.000010    0.0074      0.000010    4.55     4.55     |
 -----------------------------------------------------------------------------------------------------------
| Totals:                       3725      0.1624                                          100.00            |
 -----------------------------------------------------------------------------------------------------------

 ----------------------------------------------------------------------------------------------------------------
| libMesh Performance: Alive time=2.68763, Active time=2.627                                                     |
 ----------------------------------------------------------------------------------------------------------------
| 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.0117      0.011731    0.0167      0.016661    0.45     0.63     |
|   build_constraint_matrix()        745       0.0056      0.000008    0.0056      0.000008    0.22     0.22     |
|   build_sparsity()                 1         0.0184      0.018381    0.0459      0.045875    0.70     1.75     |
|   cnstrn_elem_mat_vec()            745       0.0183      0.000025    0.0183      0.000025    0.70     0.70     |
|   create_dof_constraints()         1         0.0590      0.059013    0.1561      0.156096    2.25     5.94     |
|   distribute_dofs()                1         0.0175      0.017502    0.0746      0.074562    0.67     2.84     |
|   dof_indices()                    19297     0.1000      0.000005    0.1000      0.000005    3.81     3.81     |
|   prepare_send_list()              1         0.0002      0.000182    0.0002      0.000182    0.01     0.01     |
|   reinit()                         1         0.0282      0.028162    0.0282      0.028162    1.07     1.07     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          2         0.0128      0.006386    0.0533      0.026644    0.49     2.03     |
|                                                                                                                |
| ExodusII_IO                                                                                                    |
|   write_nodal_data()               1         0.0581      0.058085    0.0581      0.058085    2.21     2.21     |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        1433      0.0283      0.000020    0.0283      0.000020    1.08     1.08     |
|   init_shape_functions()           689       0.0013      0.000002    0.0013      0.000002    0.05     0.05     |
|   inverse_map()                    25764     0.0597      0.000002    0.0597      0.000002    2.27     2.27     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             1433      0.0152      0.000011    0.0152      0.000011    0.58     0.58     |
|   compute_face_map()               688       0.0074      0.000011    0.0074      0.000011    0.28     0.28     |
|   init_face_shape_functions()      1         0.0000      0.000016    0.0000      0.000016    0.00     0.00     |
|   init_reference_to_physical_map() 689       0.0199      0.000029    0.0199      0.000029    0.76     0.76     |
|                                                                                                                |
| GMVIO                                                                                                          |
|   write_nodal_data()               1         0.0896      0.089639    0.0897      0.089706    3.41     3.41     |
|                                                                                                                |
| LocationMap                                                                                                    |
|   find()                           37744     0.0434      0.000001    0.0434      0.000001    1.65     1.65     |
|   init()                           1         0.0019      0.001860    0.0019      0.001860    0.07     0.07     |
|                                                                                                                |
| Mesh                                                                                                           |
|   find_neighbors()                 2         0.0673      0.033654    0.0677      0.033861    2.56     2.58     |
|   renumber_nodes_and_elem()        4         0.0051      0.001268    0.0051      0.001268    0.19     0.19     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   compute_hilbert_indices()        3         0.0632      0.021051    0.0632      0.021051    2.40     2.40     |
|   find_global_indices()            3         0.0078      0.002596    0.0771      0.025688    0.30     2.93     |
|   parallel_sort()                  3         0.0013      0.000417    0.0042      0.001414    0.05     0.16     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         2         0.0001      0.000063    0.2014      0.100687    0.00     7.67     |
|                                                                                                                |
| MeshRefinement                                                                                                 |
|   _refine_elements()               1         0.0563      0.056262    0.1761      0.176074    2.14     6.70     |
|   add_point()                      37744     0.0681      0.000002    0.1156      0.000003    2.59     4.40     |
|   make_flags_parallel_consistent() 1         0.0023      0.002288    0.0027      0.002668    0.09     0.10     |
|   make_refinement_compatible()     2         0.0019      0.000949    0.0019      0.000963    0.07     0.07     |
|                                                                                                                |
| MeshTools::Generation                                                                                          |
|   build_cube()                     1         0.0038      0.003838    0.0038      0.003838    0.15     0.15     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      2         0.4539      0.226936    0.5139      0.256929    17.28    19.56    |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      13        0.0229      0.001765    0.0242      0.001859    0.87     0.92     |
|   max(bool)                        5         0.0001      0.000025    0.0001      0.000025    0.00     0.00     |
|   max(scalar)                      201       0.0012      0.000006    0.0012      0.000006    0.05     0.05     |
|   max(vector)                      47        0.0003      0.000007    0.0009      0.000020    0.01     0.04     |
|   min(bool)                        241       0.0010      0.000004    0.0010      0.000004    0.04     0.04     |
|   min(scalar)                      194       0.2754      0.001420    0.2754      0.001420    10.48    10.48    |
|   min(vector)                      47        0.0004      0.000009    0.0013      0.000028    0.02     0.05     |
|   probe()                          63        0.0134      0.000212    0.0134      0.000212    0.51     0.51     |
|   receive()                        63        0.0003      0.000005    0.0137      0.000217    0.01     0.52     |
|   send()                           63        0.0002      0.000004    0.0002      0.000004    0.01     0.01     |
|   send_receive()                   66        0.0004      0.000006    0.0098      0.000148    0.02     0.37     |
|   sum()                            26        0.0047      0.000181    0.0052      0.000200    0.18     0.20     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           63        0.0001      0.000001    0.0001      0.000001    0.00     0.00     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         2         0.0095      0.004772    0.1727      0.086344    0.36     6.57     |
|   set_parent_processor_ids()       2         0.0025      0.001257    0.0025      0.001257    0.10     0.10     |
|                                                                                                                |
| PetscLinearSolver                                                                                              |
|   solve()                          1         0.8788      0.878796    0.8788      0.878796    33.45    33.45    |
|                                                                                                                |
| System                                                                                                         |
|   assemble()                       1         0.0882      0.088189    0.1926      0.192640    3.36     7.33     |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            128105    2.6270                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example subdomains_ex1:
*  mpirun -np 4 example-devel -d 3 -n 15 -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/subdomains/subdomains_ex1'

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

Hosted By:
SourceForge.net Logo