The source file biharmonic.h with comments:

        #ifndef __biharmonic_h__
        #define __biharmonic_h__
        
        #include "libmesh/equation_systems.h"
        #include "libmesh/serial_mesh.h"
        #include "libmesh/exodusII_io.h"
        
Bring in bits from the libMesh namespace. Just the bits we're using, since this is a header.
        using libMesh::EquationSystems;
        using libMesh::ExodusII_IO;
        using libMesh::MeshRefinement;
        using libMesh::Point;
        using libMesh::Real;
        using libMesh::UnstructuredMesh;
        
libmesh_error() and libmesh_assert() macros with a message
        #define ERROR(message)                                                                                         \
          do {                                                                                                         \
            libMesh::err << "Error: " << message << "\n";                                                              \
            libmesh_error();							                                       \
          } while(0)
        
        #define ASSERT(asserted, message) \
          do {                    \
            if(!(asserted)) {     \
              libMesh::err << "Assertion '" #asserted "' violated: " #message; \
              libmesh_error();    \
            }                     \
          } while(0)
        
        
        /**
         * The Biharmonic class encapsulates most of the data structures
         * necessary to calculate the biharmonic residual and Jacobian,
         * auxiliary quantities, to take a timestep, and to output the state --
         * biharmonic solution and vectors of auxiliary quantities.
         *
         * The main reason for this design is to have a data structure that
         * has all of the necessary data in one place, where all of the
         * calculation subroutines can access these data. Currently these data
         * are split up among several interdependent objects with no clear
         * hierarchy between them: mesh, equation system, equation system
         * bundle, residual/Jacobian calculator.
         *
         * Since no object contains all others and the data are distributed
         * among many objects, the natural control and data flow resides outside
         * of these objects and is typically implemented in main().  We,
         * however, would like to split the calculation into natural chunks --
         * subroutines -- while retaining these subroutines access to the common
         * necessary data -- biharmonic parameters, mesh and time interval
         * sizes, etc. Thus, class Biharmonic.  Finally, making Biharmonic
         * inherit from EquationSystems makes it possible to include it in the
         * most common callbacks that do not pass back a user context, but only
         * an EquationSystems object.
         */
        class Biharmonic : public EquationSystems
        {
        public:
Initial state enumeration
          enum InitialStateEnum {STRIP = 0,
        			 ROD   = 1,
        			 BALL  = 2};
        
Free energy enumeration
          enum FreeEnergyEnum {DOUBLE_WELL         = 1,
        		       DOUBLE_OBSTACLE     = 2,
        		       LOG_DOUBLE_WELL     = 3,
        		       LOG_DOUBLE_OBSTACLE = 4};
        
          /**
           * Static creation/destruction routines.  FIXME - this looks like
           * object-oriented C, can we get rid of it?
           */
          static void Create(Biharmonic** b, const Parallel::Communicator &comm);
          static void Destroy(Biharmonic** b);
        
        
          /**
           * Constructor retrieves command-line options, setting  defaults, if necessary.
           * It then builds the mesh using these options, then the equations systems around it,
           * and, finally, sets up the output.
           * We recommend that this be used through the factory Create function, which allocates
           * the mesh. In that case don't forget to call Destroy at the end, to free the mesh up.
           */
          Biharmonic(UnstructuredMesh* m);
        
        
          /**
           * Destructor
           */
          ~Biharmonic()
          {
delete _meshRefinement;
          }
        
        
Misc. getters
          bool verbose()         { return _verbose; }
          Real dt0()             { return _dt0; }
          Real dt()              { return _dt; }
        
        
Public interface functions
          void viewParameters();
          void init();
          void step(const Real& dt = -1.0);
          void output(int timestep, const Real& t, Real& o_t, bool force = false);
          void run();
        
        private:
          unsigned int  _dim, _N;
          Real _kappa, _theta, _theta_c;
          Real _tol;
          bool _growth, _degenerate, _cahn_hillard, _netforce;
          FreeEnergyEnum  _energy;
          int _log_truncation;
          bool _verbose;
          InitialStateEnum  _initialState;
          Point _initialCenter;
          Real _initialWidth;
          Real _dt0, _dt, _t0, _T, _t1;
          Real _cnWeight;


          std::string  _ofile_base, _ofile;
          ExodusII_IO* _exio;
          Real    _o_dt;
          int     _o_count;


          friend class JR;
          class JR;       // forward
          UnstructuredMesh*                       _mesh;
          MeshRefinement*                         _meshRefinement;
          JR*                                     _jr;
        };
        
        
        
        
        
        
        
        #endif // __biharmonic_h__



The source file biharmonic_jr.h with comments:

        #ifndef __biharmonic_jr_h__
        #define __biharmonic_jr_h__
        
LibMesh includes
        #include "libmesh/transient_system.h"
        #include "libmesh/nonlinear_solver.h"
        
        
Example includes
        #include "biharmonic.h"
        
Bring in bits from the libMesh namespace. Just the bits we're using, since this is a header.
        using libMesh::EquationSystems;
        using libMesh::Gradient;
        using libMesh::NonlinearImplicitSystem;
        using libMesh::Number;
        using libMesh::NumericVector;
        using libMesh::Parameters;
        using libMesh::Point;
        using libMesh::SparseMatrix;
        using libMesh::System;
        using libMesh::TransientNonlinearImplicitSystem;
        
        
        /**
         * Biharmonic's friend class definition
         */
        class Biharmonic::JR : public TransientNonlinearImplicitSystem,
        		       public NonlinearImplicitSystem::ComputeResidualandJacobian,
        		       public NonlinearImplicitSystem::ComputeBounds,
        		       public System::Initialization
        {
        public:
          /**
           * Constructor.
           */
          JR(EquationSystems& eqSys, const std::string& name, const unsigned int number);
        
          void initialize();
        
          /**
           * Static functions to be used for initialization
           */
          static Number InitialDensityBall(const Point& p, const Parameters& parameters, const std::string&, const std::string&);
          static Number InitialDensityRod(const Point& p, const Parameters& parameters, const std::string&, const std::string&);
          static Number InitialDensityStrip(const Point& p, const Parameters& parameters, const std::string&, const std::string&);
          static Gradient InitialGradientZero(const Point&, const Parameters&, const std::string&, const std::string&);
        
          /**
           * The residual and Jacobian assembly function for the Biharmonic system.
           */
          void residual_and_jacobian(const NumericVector<Number>& u,
        			     NumericVector<Number>* R,
        			     SparseMatrix<Number>* J,
        			     NonlinearImplicitSystem&);
        
        
          /**
           * Function defining the bounds of the Biharmonic system.
           */
          void bounds(NumericVector<Number>& XL,
        	      NumericVector<Number>& XU,
        	      NonlinearImplicitSystem&);
        
        private:
          Biharmonic& _biharmonic;
        };
        
        
        #endif // __biharmonic_jr_h__



The source file biharmonic.C with comments:

        #include "libmesh/mesh_generation.h"
        #include "libmesh/numeric_vector.h"
        #include "libmesh/serial_mesh.h"
        
Example includes
        #include "biharmonic.h"
        #include "biharmonic_jr.h"
        
        using namespace libMesh;
        
        void Biharmonic::Create(Biharmonic** b, const Parallel::Communicator &comm)
        {
ParallelMesh doesn't yet understand periodic BCs
          SerialMesh* mesh = new SerialMesh(comm);
          Biharmonic *biharmonic = new Biharmonic(mesh);
          *b = biharmonic;
        }
        
        
        
        
        
        void Biharmonic::Destroy(Biharmonic** b)
        {
          Biharmonic* biharmonic = *b;
          UnstructuredMesh* mesh = biharmonic->_mesh;
          delete biharmonic;
          delete mesh;
          *b = NULL;
        }
        
        
        
        void Biharmonic::viewParameters()
        {
          libMesh::out << "Biharmonic parameters:\n";
        
Print verbosity status
          if (_verbose)
            libMesh::out << "verbose mode is on\n";
          else
            libMesh::out << "verbose mode is off\n";
        
Print parameters
          libMesh::out << "mesh dimension           = " << _dim <<               "\n";
          libMesh::out << "initial linear mesh size = " << _N   <<               "\n";
          libMesh::out << "kappa                    = " << _kappa <<             "\n";
          libMesh::out << "growth                   = " << (int)_growth <<       "\n";
          libMesh::out << "degenerate               = " << (int)_degenerate <<   "\n";
          libMesh::out << "Cahn-Hillard             = " << (int)_cahn_hillard << "\n";
          libMesh::out << "netforce                 = " << (int)_netforce <<     "\n";
          libMesh::out << "energy                   = " << _energy        <<     "\n";
          libMesh::out << "tol                      = " << _tol           <<     "\n";
          libMesh::out << "theta                    = " << _theta         <<     "\n";
          libMesh::out << "theta_c                  = " << _theta_c       <<     "\n";
          libMesh::out << "log truncation           = " << _log_truncation <<    "\n";
          libMesh::out << "initial timestep size    = " << _dt0            <<    "\n";
        
          if (_initialState == STRIP)
            libMesh::out << "initial state:             strip\n";
        
          if (_initialState == ROD)
            libMesh::out << "initial state:             rod\n";
        
          if (_initialState == BALL)
            libMesh::out << "initial state:             ball\n";
        
          libMesh::out << "initial state center     = " << _initialCenter(0) << "\n";
          libMesh::out << "initial state width      = " << _initialWidth << "\n";
          libMesh::out << "initial time (min_time)  = " << _t0 << "\n";
          libMesh::out << "integration time         = " << _T  << "\n";
          libMesh::out << "final time   (max_time)  = " << _t1 << "\n";
          libMesh::out << "Crank-Nicholson weight   = " << _cnWeight << "\n";
          libMesh::out << "Output timestep          = " << _o_dt << "\n";
          libMesh::out << "Output filename base:      " <<  _ofile_base << "\n";
        }
        
        
        
        
        void Biharmonic::init()
        {
          if(_verbose)
            libMesh::out << ">>> Initializing Biharmonic\n";
        
          _dt  =  0;
          _o_count = 0;
          this->EquationSystems::init();
        
          if(_verbose)
            libMesh::out << "<<< Initializing Biharmonic\n";
        }
        
        
        
        
        
        void Biharmonic::step(const Real& dt_in)
        {
We need to update the old solution vector. The old solution vector will be the current solution vector from the previous time step. We use vector assignment. Only \p TransientSystems (and systems derived from them) contain old solutions.
          if (dt_in < 0)
            _dt = _dt0;
          else
            _dt = dt_in;
        
          *(_jr->old_local_solution) = *(_jr->current_local_solution);
        
this will localize the current solution, resulting in a current_local_solution with correct ghost values
          _jr->solve();
        }
        
        
        
        void Biharmonic::output(int timestep, const Real& t, Real& o_t, bool force)
        {
        #ifdef LIBMESH_HAVE_EXODUS_API
          if (!force && t - o_t < _o_dt)
            return;
        
          ++_o_count;
        
          if (_verbose)
            libMesh::out << "Writing state " << timestep << " at time " << t << " to file " << _ofile << "; output a total of " << _o_count << " states so far\n";
        
          _exio->write_timestep(_ofile, *this, timestep, t);
        
          if (!force)
            o_t = t;
        #endif // #ifdef LIBMESH_HAVE_EXODUS_API
        }
        
        
        
        void Biharmonic::run()
        {
          Real t = _t0, o_t = 0.0;
          int timestep = 1;
        
Force-write the initial timestep
          output(timestep,t,o_t,true);
        
          while (t < _t1)
            {
              ++timestep;
        
A pretty update message
              if (_verbose)
        	libMesh::out << "Solving for state " << timestep << ", time " << t << "\n";
        
Move biharmonic one timestep forward
              step();
        
Keep track of time and timestep
              t += _dt;
        
Output
              output(timestep,t,o_t);
            } // while(t < _t1)
        
Force-write the final timestep
          output(timestep,t,o_t,true);
        }
        
        
        
        
        
        Biharmonic::Biharmonic(UnstructuredMesh* m) :
            EquationSystems(*m),
            _mesh(m)
          {
Retrieve parameters and set defaults
            _verbose      = false; if(on_command_line("--verbose")) _verbose = true;
            _growth       = false; if(on_command_line("--growth"))       _growth = true;
            _degenerate   = false; if(on_command_line("--degenerate"))   _degenerate = true;
            _cahn_hillard = false; if(on_command_line("--cahn_hillard")) _cahn_hillard = true;
            _netforce     = false; if(on_command_line("--netforce"))     _netforce = true;
            _kappa = command_line_value("kappa", 1.0);
        
"type of energy (double well, double obstacle, logarithmic+double well, logarithmic+double obstacle)"
            std::string energy = command_line_value("energy", std::string("double_well"));
            if (energy == "double_well")
              _energy = DOUBLE_WELL;
            else if (energy == "double_obstacle")
              _energy = DOUBLE_OBSTACLE;
            else if (energy == "log_double_well")
              _energy = LOG_DOUBLE_WELL;
            else if (energy == "log_double_obstacle")
              _energy = LOG_DOUBLE_OBSTACLE;
            else
              ERROR(std::string("Unknown energy type: ") + energy);
        
            _tol     = command_line_value("tol",1.0e-8);
            _theta   = command_line_value("theta", .001);
            _theta_c = command_line_value("theta_c",1.0);
        
"order of log truncation (0=none, 2=quadratic, 3=cubic)"
            _log_truncation = command_line_value("log_truncation", 2);
        
            if (!_log_truncation)
              libMesh::out << "WARNING: no truncation is being used for the logarithmic free energy term.\nWARNING: division by zero possible!\n";
        
        
Dimension
            _dim = command_line_value("dim",1);
        
            ASSERT((_dim <= 3) && (_dim > 0), "Invalid mesh dimension");
        
Build the mesh Yes, it's better to make a coarse mesh and then refine it. We'll get to it later.
            _N = command_line_value("N", 8);
            ASSERT(_N > 0, "Invalid mesh size");
        
            switch (_dim)
              {
              case 1:
        	MeshTools::Generation::build_line(*_mesh, _N, 0.0, 1.0, EDGE2);
        	break;
              case 2:
        	MeshTools::Generation::build_square(*_mesh, _N, _N, 0.0, 1.0, 0.0, 1.0, QUAD4);
        	break;
              case 3:
        	MeshTools::Generation::build_cube(*_mesh, _N, _N, _N, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, HEX8);
        	break;
              default:
        	ASSERT((_dim <= 3) && (_dim > 0), "Invalid mesh dimension");
        	break;
              }
        
Determine the initial timestep size
            _dt0 = command_line_value("dt", 1.0/(10*_kappa*_N*_N*_N*_N));
            ASSERT(_dt0>=0, "Negative initial timestep");
        
            _t0 = command_line_value("min_time", 0.0);
            _t1 = command_line_value("max_time", _t0 + 50.0*_dt0);
            ASSERT(_t1 >= _t0, "Final time less than initial time");
            _T = _t1 - _t0;
        
            _cnWeight = command_line_value("crank_nicholson_weight", 1.0);
            ASSERT(_cnWeight <= 1 && _cnWeight >= 0, "Crank-Nicholson weight must be between 0 and 1");
        
Initial state
            _initialState = STRIP;
            std::string initialState = command_line_value("initial_state", std::string("strip"));
            if (initialState == std::string("ball"))
              _initialState = BALL;
            else if (initialState == std::string("strip"))
              _initialState = STRIP;
            else if (initialState == std::string("rod"))
              _initialState = ROD;
            else
              ERROR("Unknown initial state: neither ball nor rod nor srip");
        
            std::vector<Real> icenter;
            command_line_vector("initial_center", icenter);
        
Check that the point defining the center was in the right spatial dimension
            if (icenter.size() > _dim)
              ASSERT(icenter.size() > _dim, "Invalid dimension for the initial state center of mass");
        
Pad
            icenter.resize(3);
            for (unsigned int i = icenter.size(); i < _dim; ++i)
              icenter[i] = 0.5;
        
            for (unsigned int i = _dim; i < 3; ++i)
              icenter[i] = 0.0;
        
            _initialCenter = Point(icenter[0],icenter[1], icenter[2]);
            _initialWidth = command_line_value("initial_width", 0.125);
        
Build the main equation encapsulated in the JR (Jacobian-Residual or J(R) "jet of R") object
            _jr = &(add_system<Biharmonic::JR>(std::string("Biharmonic::JR")));
        
Output options
        #ifdef LIBMESH_HAVE_EXODUS_API
            if (on_command_line("output_base"))
              _ofile_base = command_line_value("output_base", std::string("bih"));
        
            else
              {
        	switch(_dim)
        	  {
        	  case 1:
        	    _ofile_base = std::string("bih.1");
        	    break;
        	  case 2:
        	    _ofile_base = std::string("bih.2");
        	    break;
        	  case 3:
        	    _ofile_base = std::string("bih.3");
        	    break;
        	  default:
        	    _ofile_base = std::string("bih");
        	    break;
        	  }
              }
            _ofile = _ofile_base + ".e";
            _exio = new ExodusII_IO(*_mesh);
            _o_dt = command_line_value("output_dt", 0.0);
        #endif // #ifdef LIBMESH_HAVE_EXODUS_API
          } // constructor



The source file biharmonic_jr.C with comments:

        #include "libmesh/mesh.h"
        #include "libmesh/quadrature.h"
        #include "libmesh/dense_matrix.h"
        #include "libmesh/dense_vector.h"
        #include "libmesh/sparse_matrix.h"
        #include "libmesh/fourth_error_estimators.h"
        #include "libmesh/dof_map.h"
        #include "libmesh/numeric_vector.h"
        #include "libmesh/periodic_boundaries.h"
        #include "libmesh/periodic_boundary.h"
        
Example includes
        #include "biharmonic_jr.h"
        
        using namespace libMesh;
        
        Biharmonic::JR::JR(EquationSystems& eqSys,
        		   const std::string& name_in,
        		   const unsigned int number_in) :
          TransientNonlinearImplicitSystem(eqSys,name_in,number_in),
          _biharmonic(dynamic_cast<Biharmonic&>(eqSys))
        {
Check that we can actually compute second derivatives
        #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
          ERROR("Must have second derivatives enabled");
        #endif
        
        #ifdef LIBMESH_ENABLE_PERIODIC
Add periodicity to the mesh
          DofMap& dof_map = get_dof_map();
          PeriodicBoundary xbdry(RealVectorValue(1.0, 0.0, 0.0));
        #if LIBMESH_DIM > 1
          PeriodicBoundary ybdry(RealVectorValue(0.0, 1.0, 0.0));
        #endif
        #if LIBMESH_DIM > 2
          PeriodicBoundary zbdry(RealVectorValue(0.0, 0.0, 1.0));
        #endif
        
          switch(_biharmonic._dim)
            {
            case 1:
              xbdry.myboundary = 0;
              xbdry.pairedboundary = 1;
              dof_map.add_periodic_boundary(xbdry);
              break;
        #if LIBMESH_DIM > 1
            case 2:
              xbdry.myboundary = 3;
              xbdry.pairedboundary = 1;
              dof_map.add_periodic_boundary(xbdry);
              ybdry.myboundary = 0;
              ybdry.pairedboundary = 2;
              dof_map.add_periodic_boundary(ybdry);
              break;
        #endif
        #if LIBMESH_DIM > 2
            case 3:
              xbdry.myboundary = 4;
              xbdry.pairedboundary = 2;
              dof_map.add_periodic_boundary(xbdry);
              ybdry.myboundary = 1;
              ybdry.pairedboundary = 3;
              dof_map.add_periodic_boundary(ybdry);
              zbdry.myboundary = 0;
              zbdry.pairedboundary = 5;
              dof_map.add_periodic_boundary(zbdry);
              break;
        #endif
            default:
              libmesh_error();
            }
        #endif // LIBMESH_ENABLE_PERIODIC
        
Adaptivity stuff is commented out for now... #ifndef LIBMESH_ENABLE_AMR libmesh_example_assert(false, "--enable-amr"); #else // In case we ever get around to doing mesh refinement. _biharmonic._meshRefinement = new MeshRefinement(_mesh);

// Tell the MeshRefinement object about the periodic boundaries // so that it can get heuristics like level-one conformity and unrefined // island elimination right. _biharmonic._mesh_refinement->set_periodic_boundaries_ptr(dof_map.get_periodic_boundaries()); #endif // LIBMESH_ENABLE_AMR

Adds the variable "u" to the system. u will be approximated using Hermite elements
          add_variable("u", THIRD, HERMITE);
        
Give the system an object to compute the initial state.
          attach_init_object(*this);
        
Attache the R & J calculation object
          nonlinear_solver->residual_and_jacobian_object = this;
        
Attach the bounds calculation object
          nonlinear_solver->bounds_object = this;
        }
        
        
        
        
        
        void Biharmonic::JR::initialize()
        {
          if (_biharmonic._verbose)
            libMesh::out << ">>> Initializing Biharmonic::JR\n";
        
          Parameters parameters;
          parameters.set<Point>("center") = _biharmonic._initialCenter;
          parameters.set<Real>("width")   = _biharmonic._initialWidth;
        
          if (_biharmonic._initialState == Biharmonic::BALL)
            project_solution(Biharmonic::JR::InitialDensityBall, Biharmonic::JR::InitialGradientZero, parameters);
        
          if (_biharmonic._initialState == Biharmonic::ROD)
            project_solution(Biharmonic::JR::InitialDensityRod, Biharmonic::JR::InitialGradientZero, parameters);
        
          if (_biharmonic._initialState == Biharmonic::STRIP)
            project_solution(Biharmonic::JR::InitialDensityStrip, Biharmonic::JR::InitialGradientZero, parameters);
        
both states are equal
          *(old_local_solution) = *(current_local_solution);
        
          if (_biharmonic._verbose)
            libMesh::out << "<<< Initializing Biharmonic::JR\n";
        }
        
        
        
        
        
        
        Number Biharmonic::JR::InitialDensityBall(const Point& p,
        					  const Parameters& parameters,
        					  const std::string&,
        					  const std::string&)
        {
Initialize with a ball in the middle, which is a segment in 1D, a disk in 2D and a ball in 3D.
          Point center = parameters.get<Point>("center");
          Real width = parameters.get<Real>("width");
          Point pc = p-center;
          Real r = pc.size();
          return (r < width) ? 1.0 : -0.5;
        }
        
        
        
        
        Number Biharmonic::JR::InitialDensityRod(const Point& p,
        					 const Parameters& parameters,
        					 const std::string&,
        					 const std::string&)
        {
Initialize with a rod in the middle so that we have a z-homogeneous system to model the 2D disk.
          Point center = parameters.get<Point>("center");
          Real width = parameters.get<Real>("width");
          Real r = sqrt((p(0)-center(0))*(p(0)-center(0)) + (p(1)-center(1))*(p(1)-center(1)));
          return (r < width) ? 1.0 : -0.5;
        }
        
        
        
        
        
        Number Biharmonic::JR::InitialDensityStrip(const Point& p,
        					   const Parameters& parameters,
        					   const std::string&,
        					   const std::string&)
        {
Initialize with a wide strip in the middle so that we have a yz-homogeneous system to model the 1D.
          Point center = parameters.get<Point>("center");
          Real width = parameters.get<Real>("width");
          Real r = sqrt((p(0)-center(0))*(p(0)-center(0)));
          return (r < width) ? 1.0 : -0.5;
        }
        
        
        
        
        Gradient Biharmonic::JR::InitialGradientZero(const Point&,
        					     const Parameters&,
        					     const std::string&,
        					     const std::string&)
        {
          return Gradient(0.0,0.0,0.0);
        }
        
        
        
        
        void Biharmonic::JR::residual_and_jacobian(const NumericVector<Number> &u,
        					   NumericVector<Number> *R,
        					   SparseMatrix<Number> *J,
        					   NonlinearImplicitSystem&)
        {
        #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
          if (!R && !J)
            return;
        
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 ("Biharmonic Residual and Jacobian", false);
        
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 = 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(_biharmonic._dim, fe_type));
        
Quadrature rule for numerical integration. With 2D triangles, the Clough quadrature rule puts a Gaussian quadrature rule on each of the 3 subelements
          AutoPtr<QBase> qrule(fe_type.default_quadrature_rule(_biharmonic._dim));
        
Tell the finite element object to use our quadrature rule.
          fe->attach_quadrature_rule (qrule.get());
        
Here we define some references to element-specific data that will be used to assemble the linear system. We begin with the element Jacobian * quadrature weight at each integration point.
          const std::vector<Real>& JxW = fe->get_JxW();
        
The element shape functions evaluated at the quadrature points.
          const std::vector<std::vector<Real> >& phi = fe->get_phi();
        
The element shape functions' derivatives evaluated at the quadrature points.
          const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
        
The element shape functions' second derivatives evaluated at the quadrature points.
          const std::vector<std::vector<RealTensor> >& d2phi = fe->get_d2phi();
        
For efficiency we will compute shape function laplacians n times, not n^2
          std::vector<Real> Laplacian_phi_qp;
        
Define data structures to contain the element matrix and right-hand-side vector contribution. Following basic finite element terminology we will denote these "Je" and "Re". More detail is in example 3.
          DenseMatrix<Number> Je;
          DenseVector<Number> Re;
        
This vector will hold the degree of freedom indices for the element. These define where in the global system the element degrees of freedom get mapped.
          std::vector<dof_id_type> dof_indices;
        
Old solution
          const NumericVector<Number>& u_old = *old_local_solution;
        
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     = _biharmonic._mesh->active_local_elements_begin();
          const MeshBase::const_element_iterator end_el = _biharmonic._mesh->active_local_elements_end();
        
          for ( ; el != end_el; ++el) {
Store a pointer to the element we are currently working on. This allows for nicer syntax later.
            const Elem* elem = *el;
        
Get the degree of freedom indices for the current element. These define where in the global matrix and right-hand-side this element will contribute to.
            dof_map.dof_indices (elem, dof_indices);
        
Compute the element-specific data for the current element. This involves computing the location of the quadrature points (q_point) and the shape function values/derivatives (phi, dphi,d2phi) for the current element.
            fe->reinit (elem);
        
Zero the element matrix, the right-hand side and the Laplacian matrix before summing them.
            if (J)
              Je.resize(dof_indices.size(), dof_indices.size());
        
            if (R)
              Re.resize(dof_indices.size());
        
            Laplacian_phi_qp.resize(dof_indices.size());
        
            for (unsigned int qp=0; qp<qrule->n_points(); qp++)
              {
AUXILIARY QUANTITIES: Residual and Jacobian share a few calculations: at the very least, in the case of interfacial energy only with a constant mobility, both calculations use Laplacian_phi_qp; more is shared the case of a concentration-dependent mobility and bulk potentials.
                Number u_qp = 0.0, u_old_qp = 0.0, Laplacian_u_qp = 0.0, Laplacian_u_old_qp = 0.0;
        	Gradient grad_u_qp(0.0,0.0,0.0), grad_u_old_qp(0.0,0.0,0.0);
        	Number M_qp = 1.0, M_old_qp = 1.0, M_prime_qp = 0.0, M_prime_old_qp = 0.0;
        
        	for (unsigned int i=0; i<phi.size(); i++)
        	  {
        	    Laplacian_phi_qp[i] = d2phi[i][qp](0,0);
        	    grad_u_qp(0) += u(dof_indices[i])*dphi[i][qp](0);
        	    grad_u_old_qp(0) += u_old(dof_indices[i])*dphi[i][qp](0);
        
        	    if (_biharmonic._dim > 1)
        	      {
        		Laplacian_phi_qp[i] += d2phi[i][qp](1,1);
        		grad_u_qp(1) += u(dof_indices[i])*dphi[i][qp](1);
        		grad_u_old_qp(1) += u_old(dof_indices[i])*dphi[i][qp](1);
        	      }
        	    if (_biharmonic._dim > 2)
        	      {
        		Laplacian_phi_qp[i] += d2phi[i][qp](2,2);
        		grad_u_qp(2) += u(dof_indices[i])*dphi[i][qp](2);
        		grad_u_old_qp(2) += u_old(dof_indices[i])*dphi[i][qp](2);
        	      }
        	    u_qp     += phi[i][qp]*u(dof_indices[i]);
        	    u_old_qp += phi[i][qp]*u_old(dof_indices[i]);
        	    Laplacian_u_qp     += Laplacian_phi_qp[i]*u(dof_indices[i]);
        	    Laplacian_u_old_qp += Laplacian_phi_qp[i]*u_old(dof_indices[i]);
        	  } // for i
        
        	if (_biharmonic._degenerate)
        	  {
        	    M_qp           = 1.0 - u_qp*u_qp;
        	    M_old_qp       = 1.0 - u_old_qp*u_old_qp;
        	    M_prime_qp     = -2.0*u_qp;
        	    M_prime_old_qp = -2.0*u_old_qp;
        	  }
        
ELEMENT RESIDUAL AND JACOBIAN
                for (unsigned int i=0; i<phi.size(); i++)
        	  {
RESIDUAL
                    if (R)
        	      {
        		Number ri = 0.0, ri_old = 0.0;
        		ri     -= Laplacian_phi_qp[i]*M_qp*_biharmonic._kappa*Laplacian_u_qp;
        		ri_old -= Laplacian_phi_qp[i]*M_old_qp*_biharmonic._kappa*Laplacian_u_old_qp;
        
        		if (_biharmonic._degenerate)
        		  {
        		    ri       -= (dphi[i][qp]*grad_u_qp)*M_prime_qp*(_biharmonic._kappa*Laplacian_u_qp);
        		    ri_old   -= (dphi[i][qp]*grad_u_old_qp)*M_prime_old_qp*(_biharmonic._kappa*Laplacian_u_old_qp);
        		  }
        
        		if (_biharmonic._cahn_hillard)
        		  {
        		    if (_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
        		      {
        			ri += Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp;
        			ri_old += Laplacian_phi_qp[i]*M_old_qp*_biharmonic._theta_c*(u_old_qp*u_old_qp - 1.0)*u_old_qp;
        			if (_biharmonic._degenerate)
        			  {
        			    ri     += (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp;
        			    ri_old += (dphi[i][qp]*grad_u_old_qp)*M_prime_old_qp*_biharmonic._theta_c*(u_old_qp*u_old_qp - 1.0)*u_old_qp;
        			  }
        		      }// if(_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
        
        		    if (_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
        		      {
        			ri -= Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*u_qp;
        			ri_old -= Laplacian_phi_qp[i]*M_old_qp*_biharmonic._theta_c*u_old_qp;
        			if (_biharmonic._degenerate)
        			  {
        			    ri     -= (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*u_qp;
        			    ri_old -= (dphi[i][qp]*grad_u_old_qp)*M_prime_old_qp*_biharmonic._theta_c*u_old_qp;
        			  }
        		      } // if(_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
        
        		    if (_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
        		      {
        			switch(_biharmonic._log_truncation)
        			  {
        			  case 2:
        			    break;
        			  case 3:
        			    break;
        			  default:
        			    break;
        			  }// switch(_biharmonic._log_truncation)
        		      }// if(_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
        		  }// if(_biharmonic._cahn_hillard)
        		Re(i) += JxW[qp]*((u_qp-u_old_qp)*phi[i][qp]-_biharmonic._dt*0.5*((2.0-_biharmonic._cnWeight)*ri + _biharmonic._cnWeight*ri_old));
        	      } // if (R)
        
JACOBIAN
                    if (J)
        	      {
        		Number M_prime_prime_qp = 0.0;
        		if(_biharmonic._degenerate) M_prime_prime_qp = -2.0;
        		for (unsigned int j=0; j<phi.size(); j++)
        		  {
        		    Number ri_j = 0.0;
        		    ri_j -= Laplacian_phi_qp[i]*M_qp*_biharmonic._kappa*Laplacian_phi_qp[j];
        		    if (_biharmonic._degenerate)
        		      {
        			ri_j -=
        			  Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._kappa*Laplacian_u_qp               +
        			  (dphi[i][qp]*dphi[j][qp])*M_prime_qp*(_biharmonic._kappa*Laplacian_u_qp)                  +
        			  (dphi[i][qp]*grad_u_qp)*(M_prime_prime_qp*phi[j][qp])*(_biharmonic._kappa*Laplacian_u_qp) +
        			  (dphi[i][qp]*grad_u_qp)*(M_prime_qp)*(_biharmonic._kappa*Laplacian_phi_qp[j]);
        		      }
        
        		    if (_biharmonic._cahn_hillard)
        		      {
        			if(_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
        			  {
        			    ri_j +=
        			      Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp +
        			      Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*(3.0*u_qp*u_qp - 1.0)*phi[j][qp]        +
        			      (dphi[i][qp]*dphi[j][qp])*M_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp      +
        			      (dphi[i][qp]*grad_u_qp)*M_prime_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp  +
        			      (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*(3.0*u_qp*u_qp - 1.0)*phi[j][qp];
        			  }// if(_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
        
        			if (_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
        			  {
        			    ri_j -=
        			      Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._theta_c*u_qp                   +
        			      Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*phi[j][qp]                              +
        			      (dphi[i][qp]*dphi[j][qp])*M_prime_qp*_biharmonic._theta_c*u_qp                        +
        			      (dphi[i][qp]*grad_u_qp)*M_prime_prime_qp*_biharmonic._theta_c*u_qp                    +
        			      (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*phi[j][qp];
        			  } // if(_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
        
        			if (_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
        			  {
        			    switch(_biharmonic._log_truncation)
        			      {
        			      case 2:
        				break;
        			      case 3:
        				break;
        			      default:
        				break;
        			      }// switch(_biharmonic._log_truncation)
        			  }// if(_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
        		      }// if(_biharmonic._cahn_hillard)
        		    Je(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] - 0.5*_biharmonic._dt*(2.0-_biharmonic._cnWeight)*ri_j);
        		  } // for j
        	      } // if (J)
        	  } // for i
              } // for qp
        
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
            if (R)
              {
If the mesh has hanging nodes (e.g., as a result of refinement), those need to be constrained.
                dof_map.constrain_element_vector(Re, dof_indices);
        	R->add_vector(Re, dof_indices);
              }
        
            if (J)
              {
If the mesh has hanging nodes (e.g., as a result of refinement), those need to be constrained.
                dof_map.constrain_element_matrix(Je, dof_indices);
        	J->add_matrix(Je, dof_indices);
              }
          } // for el
        #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
        }
        
        
        
        
        
        void Biharmonic::JR::bounds(NumericVector<Number> &XL, NumericVector<Number>& XU, NonlinearImplicitSystem&)
        {
sys is actually ignored, since it should be the same as *this.

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 ("Biharmonic bounds", false);
        
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 = 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(_biharmonic._dim, fe_type));
        
Define data structures to contain the bound vectors contributions.
          DenseVector<Number> XLe, XUe;
        
These 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;
        
          MeshBase::const_element_iterator       el     = _biharmonic._mesh->active_local_elements_begin();
          const MeshBase::const_element_iterator end_el = _biharmonic._mesh->active_local_elements_end();
        
          for ( ; el != end_el; ++el)
            {
Extract the shape function to be evaluated at the nodes
              const std::vector<std::vector<Real> >& phi = fe->get_phi();
        
Get the degree of freedom indices for the current element. They are in 1-1 correspondence with shape functions phi and define where in the global vector this element will.
              dof_map.dof_indices (*el, dof_indices);
        
Resize the local bounds vectors (zeroing them out in the process).
              XLe.resize(dof_indices.size());
              XUe.resize(dof_indices.size());
        
Extract the element node coordinates in the reference frame
              std::vector<Point> nodes;
              fe->get_refspace_nodes((*el)->type(), nodes);
        
Evaluate the shape functions at the nodes
              fe->reinit(*el, &nodes);
        
Construct the bounds based on the value of the i-th phi at the nodes. Observe that this doesn't really work in general: we rely on the fact that for Hermite elements each shape function is nonzero at most at a single node. More generally the bounds must be constructed by inspecting a "mass-like" matrix (m_{ij}) of the shape functions (i) evaluated at their corresponding nodes (j). The constraints imposed on the dofs (d_i) are then are -1 \leq \sum_i d_i m_{ij} \leq 1, since \sum_i d_i m_{ij} is the value of the solution at the j-th node. Auxiliary variables will need to be introduced to reduce this to a "box" constraint. Additional complications will arise since m might be singular (as is the case for Hermite, which, however, is easily handled by inspection).
              for (unsigned int i=0; i<phi.size(); ++i)
        	{
FIXME: should be able to define INF and pass it to the solve
                  Real infinity = 1.0e20;
        	  Real bound = infinity;
        	  for(unsigned int j = 0; j < nodes.size(); ++j) {
        	    if(phi[i][j]) {
        	      bound = 1.0/fabs(phi[i][j]);
        	      break;
        	    }
        	  }
        
The value of the solution at this node must be between 1.0 and -1.0. Based on the value of phi(i)(i) the nodal coordinate must be between 1.0/phi(i)(i) and its negative.
                  XLe(i) = -bound;
        	  XUe(i) = bound;
        	}
The element bound vectors are now built for this element. Insert them into the global vectors, potentially overwriting the same dof contributions from other elements: no matter -- the bounds are always -1.0 and 1.0.
              XL.insert(XLe, dof_indices);
              XU.insert(XUe, dof_indices);
            }
        }



The source file miscellaneous_ex7.C with comments:

        #include "biharmonic.h"
        
        
Bring in everything from the libMesh namespace
        using namespace libMesh;
        
Print usage information if requested on command line
        void print_help(int argc, char** argv);
        
        int main(int argc, char** argv)
        {
Initialize libMesh.
          LibMeshInit init (argc, argv);
          if (on_command_line("--help"))
            print_help(argc, argv);
          else
            {
        #if !defined(LIBMESH_ENABLE_SECOND_DERIVATIVES)
              libmesh_example_assert(false, "--enable-second");
        #elif !defined(LIBMESH_ENABLE_PERIODIC)
              libmesh_example_assert(false, "--enable-periodic");
        #endif
        
This is a PETSc-specific solver
              libmesh_example_assert(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
        
              const int dim = command_line_value("dim",1);
        
Skip higher-dimensional examples on a lower-dimensional libMesh build
              libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");
        
              Biharmonic* biharmonic;
              Biharmonic::Create(&biharmonic,init.comm());
              biharmonic->viewParameters();
              biharmonic->init();
              biharmonic->run();
              Biharmonic::Destroy(&biharmonic);
            }
          return 0;
        }
        
        
        
        
        
        void print_help(int, char** argv)
        {
          libMesh::out << "This example solves the Cahn-Hillard equation with chemical potential f:\n"
        	       << "    u_t = \\div(M(u)\\grad f(u))\n"
        	       << "Here we have\n"
        	       << "    u, -1 <= u <= 1        -- relative concentration (difference of two concentrations in a binary mixture) \n"
        	       << "    M, M >= 0              -- mobility of the mixture\n"
        	       << "    f = \\delta E/\\delta u  -- variational derivative of the free energy functional E\n"
        	       << "    E = \\int[\\kappa/2 |\\grac u|^ + g(u)]\n"
        	       << "where the gradient term is the interfacial energy density with \\kappa quantifying the energy of the interface,\n"
        	       << "and g(u) is the bulk energy density\n"
        	       << "    g(u) = \\theta L(u) + \\theta_c W(u),\n"
        	       << "L(u) is the (optional, in this model) logarithmic term corresponding to the entropy of the mixture:\n"
        	       << "    L(u) = (\\theta/2)[(1+u)\\ln((1+u)/2) + (1-u)\\ln((1-u)/2)],\n"
        	       << "where \\theta is related to the Boltzmann factor k_B T - a proxy for the absolute temperature T.\n"
        	       << "L can be optionally approximated ('truncated') using a quadratic or a cubic polynomial on [-1,1]\n"
        	       << "W(u) is the (optional, in this model) potential promoting demixing.  It can take the form of \n"
        	       << "a 'double well' potential\n"
        	       << "    W(u) = \\theta_c (u^4/4 - u^2/2),\n"
        	       << "         or \n"
        	       << "a 'double obstacle' potential\n"
        	       << "    W(u) = (\\theta_c/2)(1-u^2),\n"
        	       << "where \\theta_c is the critical 'temperature'.\n"
        	       << "Finally, mobility M can be constant of 'degenerate', by which we mean that M is varying with u and \n"
        	       << "vanishing (degenerating) whenever u reaches critical values +1 or -1:\n"
        	       << "    M(u) = 1.0\n"
        	       << "      or\n"
        	       << "    M(u) = (1.0 - u^2)\n"
        	       << "Degenerate mobility should generally be used only in conjunction with logarithmic free energy terms.\n\n"
        	       << "The equation is solved on a periodic domain (in 1D, 2D or 3D)\n"
        	       << "using a Galerkin formulation with C^1 elements approximating the H^2_{per} function space.\n\n"
        	       << "\n-----------\n"
        	       << "COMPILING: "
        	       << "\n-----------\n"
        	       << "Compile as follows (assuming libmesh has been built): \n"
        	       << "METHOD=<method> make \n"
        	       << "where <method> is dbg or opt.\n"
        	       << "\n-----------\n"
        	       << "HELP:        "
        	       << "\n-----------\n"
        	       << "Print this help message:\n"
        	       << argv[0] << " --help\n"
        	       << "\n-----------\n"
        	       << "RUNNING:     "
        	       << "\n-----------\n"
        	       << "Run in serial with build METHOD <method> as follows:\n"
        	       << "\n"
        	       << argv[0] << "\n"
        	       << "               [--verbose] dim=<1|2|3> N=<number_of_linear_elements> \n"
        	       << "               kappa=<kappa_value> growth=<yes|no> degenerate=<yes|no> [--cahn-hillard]                                           \n"
        	       << "               [--netforce]  energy=<double_well|double_obstacle|log_double_well|log_double_obstacle>  log_truncation_order=<2|3> \n"
        	       << "               theta=<theta_value> theta_c=<theta_c_value>                                                                        \n"
        	       << "               initial_state=<ball|rod|strip> initial_center='x [y [z]]' initial_width=<width>                                    \n"
        	       << "               min_time=<initial_time> max_time=<final_time> dt=<timestep_size> crank_nicholson_weight=<between_0_and_1>          \n"
        	       << "               output_base=<base_filename> output_dt=<output_timestep_size> [--use-petsc-dm -snes_type virs]                      \n"
        	       << "\n"
        	       << argv[0] << " --verbose \n"
        	       << "is a pretty good start.\n"
        	       << "\nModeling a 1D system with 2D or 3D (for a strip the second and third components of the center are immaterial):\n"
        	       << argv[0]<< " --verbose dim=1 N=1024 initial_state=strip initial_center=0.5 initial_width=0.1 dt=1e-10 max_time=1e-6\n"
        	       << argv[0]<< " --verbose dim=2 N=64   initial_state=strip initial_center=0.5 initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
        	       << argv[0]<< " --verbose dim=3 N=32   initial_state=strip initial_center=0.5 initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
        	       << "\n"
        	       << "Modeling a 2D system with 3D (for a rod the third component of the center is immaterial) \n"
        	       << argv[0]<< " --verbose dim=2 N=64   initial_state=rod initial_center='0.5 0.5' initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
        	       << argv[0]<< " --verbose dim=3 N=32   initial_state=rod initial_center='0.5 0.5' initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
        	       << "\n"
        	       << "A 3D system with an initial ball in the center\n"
        	       << argv[0] << " --verbose dim=3 N=32   initial_state=ball initial_center='0.5 0.5 0.5' initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
        	       << "\n"
        	       << "Add --use-petsc-dm -snes_type virs to run the variational inequality version that ensures the solution is between -1.0 and 1.0 at all times.\n\n"
        	       << std::endl;
        }



The source file biharmonic.h without comments:

 
  #ifndef __biharmonic_h__
  #define __biharmonic_h__
  
  #include "libmesh/equation_systems.h"
  #include "libmesh/serial_mesh.h"
  #include "libmesh/exodusII_io.h"
  
  using libMesh::EquationSystems;
  using libMesh::ExodusII_IO;
  using libMesh::MeshRefinement;
  using libMesh::Point;
  using libMesh::Real;
  using libMesh::UnstructuredMesh;
  
  #define ERROR(message)                                                                                         \
    do {                                                                                                         \
      libMesh::err << "Error: " << message << "\n";                                                              \
      libmesh_error();							                                       \
    } while(0)
  
  #define ASSERT(asserted, message) \
    do {                    \
      if(!(asserted)) {     \
        libMesh::err << "Assertion '" #asserted "' violated: " #message; \
        libmesh_error();    \
      }                     \
    } while(0)
  
  
  /**
   * The Biharmonic class encapsulates most of the data structures
   * necessary to calculate the biharmonic residual and Jacobian,
   * auxiliary quantities, to take a timestep, and to output the state --
   * biharmonic solution and vectors of auxiliary quantities.
   *
   * The main reason for this design is to have a data structure that
   * has all of the necessary data in one place, where all of the
   * calculation subroutines can access these data. Currently these data
   * are split up among several interdependent objects with no clear
   * hierarchy between them: mesh, equation system, equation system
   * bundle, residual/Jacobian calculator.
   *
   * Since no object contains all others and the data are distributed
   * among many objects, the natural control and data flow resides outside
   * of these objects and is typically implemented in main().  We,
   * however, would like to split the calculation into natural chunks --
   * subroutines -- while retaining these subroutines access to the common
   * necessary data -- biharmonic parameters, mesh and time interval
   * sizes, etc. Thus, class Biharmonic.  Finally, making Biharmonic
   * inherit from EquationSystems makes it possible to include it in the
   * most common callbacks that do not pass back a user context, but only
   * an EquationSystems object.
   */
  class Biharmonic : public EquationSystems
  {
  public:
    enum InitialStateEnum {STRIP = 0,
  			 ROD   = 1,
  			 BALL  = 2};
  
    enum FreeEnergyEnum {DOUBLE_WELL         = 1,
  		       DOUBLE_OBSTACLE     = 2,
  		       LOG_DOUBLE_WELL     = 3,
  		       LOG_DOUBLE_OBSTACLE = 4};
  
    /**
     * Static creation/destruction routines.  FIXME - this looks like
     * object-oriented C, can we get rid of it?
     */
    static void Create(Biharmonic** b, const Parallel::Communicator &comm);
    static void Destroy(Biharmonic** b);
  
  
    /**
     * Constructor retrieves command-line options, setting  defaults, if necessary.
     * It then builds the mesh using these options, then the equations systems around it,
     * and, finally, sets up the output.
     * We recommend that this be used through the factory Create function, which allocates
     * the mesh. In that case don't forget to call Destroy at the end, to free the mesh up.
     */
    Biharmonic(UnstructuredMesh* m);
  
  
    /**
     * Destructor
     */
    ~Biharmonic()
    {
    }
  
  
    bool verbose()         { return _verbose; }
    Real dt0()             { return _dt0; }
    Real dt()              { return _dt; }
  
  
    void viewParameters();
    void init();
    void step(const Real& dt = -1.0);
    void output(int timestep, const Real& t, Real& o_t, bool force = false);
    void run();
  
  private:
    unsigned int  _dim, _N;
    Real _kappa, _theta, _theta_c;
    Real _tol;
    bool _growth, _degenerate, _cahn_hillard, _netforce;
    FreeEnergyEnum  _energy;
    int _log_truncation;
    bool _verbose;
    InitialStateEnum  _initialState;
    Point _initialCenter;
    Real _initialWidth;
    Real _dt0, _dt, _t0, _T, _t1;
    Real _cnWeight;
    std::string  _ofile_base, _ofile;
    ExodusII_IO* _exio;
    Real    _o_dt;
    int     _o_count;
    friend class JR;
    class JR;       // forward
    UnstructuredMesh*                       _mesh;
    MeshRefinement*                         _meshRefinement;
    JR*                                     _jr;
  };
  
  
  
  
  
  
  
  #endif // __biharmonic_h__



The source file biharmonic_jr.h without comments:

 
  #ifndef __biharmonic_jr_h__
  #define __biharmonic_jr_h__
  
  #include "libmesh/transient_system.h"
  #include "libmesh/nonlinear_solver.h"
  
  
  #include "biharmonic.h"
  
  using libMesh::EquationSystems;
  using libMesh::Gradient;
  using libMesh::NonlinearImplicitSystem;
  using libMesh::Number;
  using libMesh::NumericVector;
  using libMesh::Parameters;
  using libMesh::Point;
  using libMesh::SparseMatrix;
  using libMesh::System;
  using libMesh::TransientNonlinearImplicitSystem;
  
  
  /**
   * Biharmonic's friend class definition
   */
  class Biharmonic::JR : public TransientNonlinearImplicitSystem,
  		       public NonlinearImplicitSystem::ComputeResidualandJacobian,
  		       public NonlinearImplicitSystem::ComputeBounds,
  		       public System::Initialization
  {
  public:
    /**
     * Constructor.
     */
    JR(EquationSystems& eqSys, const std::string& name, const unsigned int number);
  
    void initialize();
  
    /**
     * Static functions to be used for initialization
     */
    static Number InitialDensityBall(const Point& p, const Parameters& parameters, const std::string&, const std::string&);
    static Number InitialDensityRod(const Point& p, const Parameters& parameters, const std::string&, const std::string&);
    static Number InitialDensityStrip(const Point& p, const Parameters& parameters, const std::string&, const std::string&);
    static Gradient InitialGradientZero(const Point&, const Parameters&, const std::string&, const std::string&);
  
    /**
     * The residual and Jacobian assembly function for the Biharmonic system.
     */
    void residual_and_jacobian(const NumericVector<Number>& u,
  			     NumericVector<Number>* R,
  			     SparseMatrix<Number>* J,
  			     NonlinearImplicitSystem&);
  
  
    /**
     * Function defining the bounds of the Biharmonic system.
     */
    void bounds(NumericVector<Number>& XL,
  	      NumericVector<Number>& XU,
  	      NonlinearImplicitSystem&);
  
  private:
    Biharmonic& _biharmonic;
  };
  
  
  #endif // __biharmonic_jr_h__



The source file biharmonic.C without comments:

 
  #include "libmesh/mesh_generation.h"
  #include "libmesh/numeric_vector.h"
  #include "libmesh/serial_mesh.h"
  
  #include "biharmonic.h"
  #include "biharmonic_jr.h"
  
  using namespace libMesh;
  
  void Biharmonic::Create(Biharmonic** b, const Parallel::Communicator &comm)
  {
    SerialMesh* mesh = new SerialMesh(comm);
    Biharmonic *biharmonic = new Biharmonic(mesh);
    *b = biharmonic;
  }
  
  
  
  
  
  void Biharmonic::Destroy(Biharmonic** b)
  {
    Biharmonic* biharmonic = *b;
    UnstructuredMesh* mesh = biharmonic->_mesh;
    delete biharmonic;
    delete mesh;
    *b = NULL;
  }
  
  
  
  void Biharmonic::viewParameters()
  {
    libMesh::out << "Biharmonic parameters:\n";
  
    if (_verbose)
      libMesh::out << "verbose mode is on\n";
    else
      libMesh::out << "verbose mode is off\n";
  
    libMesh::out << "mesh dimension           = " << _dim <<               "\n";
    libMesh::out << "initial linear mesh size = " << _N   <<               "\n";
    libMesh::out << "kappa                    = " << _kappa <<             "\n";
    libMesh::out << "growth                   = " << (int)_growth <<       "\n";
    libMesh::out << "degenerate               = " << (int)_degenerate <<   "\n";
    libMesh::out << "Cahn-Hillard             = " << (int)_cahn_hillard << "\n";
    libMesh::out << "netforce                 = " << (int)_netforce <<     "\n";
    libMesh::out << "energy                   = " << _energy        <<     "\n";
    libMesh::out << "tol                      = " << _tol           <<     "\n";
    libMesh::out << "theta                    = " << _theta         <<     "\n";
    libMesh::out << "theta_c                  = " << _theta_c       <<     "\n";
    libMesh::out << "log truncation           = " << _log_truncation <<    "\n";
    libMesh::out << "initial timestep size    = " << _dt0            <<    "\n";
  
    if (_initialState == STRIP)
      libMesh::out << "initial state:             strip\n";
  
    if (_initialState == ROD)
      libMesh::out << "initial state:             rod\n";
  
    if (_initialState == BALL)
      libMesh::out << "initial state:             ball\n";
  
    libMesh::out << "initial state center     = " << _initialCenter(0) << "\n";
    libMesh::out << "initial state width      = " << _initialWidth << "\n";
    libMesh::out << "initial time (min_time)  = " << _t0 << "\n";
    libMesh::out << "integration time         = " << _T  << "\n";
    libMesh::out << "final time   (max_time)  = " << _t1 << "\n";
    libMesh::out << "Crank-Nicholson weight   = " << _cnWeight << "\n";
    libMesh::out << "Output timestep          = " << _o_dt << "\n";
    libMesh::out << "Output filename base:      " <<  _ofile_base << "\n";
  }
  
  
  
  
  void Biharmonic::init()
  {
    if(_verbose)
      libMesh::out << ">>> Initializing Biharmonic\n";
  
    _dt  =  0;
    _o_count = 0;
    this->EquationSystems::init();
  
    if(_verbose)
      libMesh::out << "<<< Initializing Biharmonic\n";
  }
  
  
  
  
  
  void Biharmonic::step(const Real& dt_in)
  {
    if (dt_in < 0)
      _dt = _dt0;
    else
      _dt = dt_in;
  
    *(_jr->old_local_solution) = *(_jr->current_local_solution);
  
    _jr->solve();
  }
  
  
  
  void Biharmonic::output(int timestep, const Real& t, Real& o_t, bool force)
  {
  #ifdef LIBMESH_HAVE_EXODUS_API
    if (!force && t - o_t < _o_dt)
      return;
  
    ++_o_count;
  
    if (_verbose)
      libMesh::out << "Writing state " << timestep << " at time " << t << " to file " << _ofile << "; output a total of " << _o_count << " states so far\n";
  
    _exio->write_timestep(_ofile, *this, timestep, t);
  
    if (!force)
      o_t = t;
  #endif // #ifdef LIBMESH_HAVE_EXODUS_API
  }
  
  
  
  void Biharmonic::run()
  {
    Real t = _t0, o_t = 0.0;
    int timestep = 1;
  
    output(timestep,t,o_t,true);
  
    while (t < _t1)
      {
        ++timestep;
  
        if (_verbose)
  	libMesh::out << "Solving for state " << timestep << ", time " << t << "\n";
  
        step();
  
        t += _dt;
  
        output(timestep,t,o_t);
      } // while(t < _t1)
  
    output(timestep,t,o_t,true);
  }
  
  
  
  
  
  Biharmonic::Biharmonic(UnstructuredMesh* m) :
      EquationSystems(*m),
      _mesh(m)
    {
      _verbose      = false; if(on_command_line("--verbose")) _verbose = true;
      _growth       = false; if(on_command_line("--growth"))       _growth = true;
      _degenerate   = false; if(on_command_line("--degenerate"))   _degenerate = true;
      _cahn_hillard = false; if(on_command_line("--cahn_hillard")) _cahn_hillard = true;
      _netforce     = false; if(on_command_line("--netforce"))     _netforce = true;
      _kappa = command_line_value("kappa", 1.0);
  
      std::string energy = command_line_value("energy", std::string("double_well"));
      if (energy == "double_well")
        _energy = DOUBLE_WELL;
      else if (energy == "double_obstacle")
        _energy = DOUBLE_OBSTACLE;
      else if (energy == "log_double_well")
        _energy = LOG_DOUBLE_WELL;
      else if (energy == "log_double_obstacle")
        _energy = LOG_DOUBLE_OBSTACLE;
      else
        ERROR(std::string("Unknown energy type: ") + energy);
  
      _tol     = command_line_value("tol",1.0e-8);
      _theta   = command_line_value("theta", .001);
      _theta_c = command_line_value("theta_c",1.0);
  
      _log_truncation = command_line_value("log_truncation", 2);
  
      if (!_log_truncation)
        libMesh::out << "WARNING: no truncation is being used for the logarithmic free energy term.\nWARNING: division by zero possible!\n";
  
  
      _dim = command_line_value("dim",1);
  
      ASSERT((_dim <= 3) && (_dim > 0), "Invalid mesh dimension");
  
      _N = command_line_value("N", 8);
      ASSERT(_N > 0, "Invalid mesh size");
  
      switch (_dim)
        {
        case 1:
  	MeshTools::Generation::build_line(*_mesh, _N, 0.0, 1.0, EDGE2);
  	break;
        case 2:
  	MeshTools::Generation::build_square(*_mesh, _N, _N, 0.0, 1.0, 0.0, 1.0, QUAD4);
  	break;
        case 3:
  	MeshTools::Generation::build_cube(*_mesh, _N, _N, _N, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, HEX8);
  	break;
        default:
  	ASSERT((_dim <= 3) && (_dim > 0), "Invalid mesh dimension");
  	break;
        }
  
      _dt0 = command_line_value("dt", 1.0/(10*_kappa*_N*_N*_N*_N));
      ASSERT(_dt0>=0, "Negative initial timestep");
  
      _t0 = command_line_value("min_time", 0.0);
      _t1 = command_line_value("max_time", _t0 + 50.0*_dt0);
      ASSERT(_t1 >= _t0, "Final time less than initial time");
      _T = _t1 - _t0;
  
      _cnWeight = command_line_value("crank_nicholson_weight", 1.0);
      ASSERT(_cnWeight <= 1 && _cnWeight >= 0, "Crank-Nicholson weight must be between 0 and 1");
  
      _initialState = STRIP;
      std::string initialState = command_line_value("initial_state", std::string("strip"));
      if (initialState == std::string("ball"))
        _initialState = BALL;
      else if (initialState == std::string("strip"))
        _initialState = STRIP;
      else if (initialState == std::string("rod"))
        _initialState = ROD;
      else
        ERROR("Unknown initial state: neither ball nor rod nor srip");
  
      std::vector<Real> icenter;
      command_line_vector("initial_center", icenter);
  
      if (icenter.size() > _dim)
        ASSERT(icenter.size() > _dim, "Invalid dimension for the initial state center of mass");
  
      icenter.resize(3);
      for (unsigned int i = icenter.size(); i < _dim; ++i)
        icenter[i] = 0.5;
  
      for (unsigned int i = _dim; i < 3; ++i)
        icenter[i] = 0.0;
  
      _initialCenter = Point(icenter[0],icenter[1], icenter[2]);
      _initialWidth = command_line_value("initial_width", 0.125);
  
      _jr = &(add_system<Biharmonic::JR>(std::string("Biharmonic::JR")));
  
  #ifdef LIBMESH_HAVE_EXODUS_API
      if (on_command_line("output_base"))
        _ofile_base = command_line_value("output_base", std::string("bih"));
  
      else
        {
  	switch(_dim)
  	  {
  	  case 1:
  	    _ofile_base = std::string("bih.1");
  	    break;
  	  case 2:
  	    _ofile_base = std::string("bih.2");
  	    break;
  	  case 3:
  	    _ofile_base = std::string("bih.3");
  	    break;
  	  default:
  	    _ofile_base = std::string("bih");
  	    break;
  	  }
        }
      _ofile = _ofile_base + ".e";
      _exio = new ExodusII_IO(*_mesh);
      _o_dt = command_line_value("output_dt", 0.0);
  #endif // #ifdef LIBMESH_HAVE_EXODUS_API
    } // constructor



The source file biharmonic_jr.C without comments:

 
  #include "libmesh/mesh.h"
  #include "libmesh/quadrature.h"
  #include "libmesh/dense_matrix.h"
  #include "libmesh/dense_vector.h"
  #include "libmesh/sparse_matrix.h"
  #include "libmesh/fourth_error_estimators.h"
  #include "libmesh/dof_map.h"
  #include "libmesh/numeric_vector.h"
  #include "libmesh/periodic_boundaries.h"
  #include "libmesh/periodic_boundary.h"
  
  #include "biharmonic_jr.h"
  
  using namespace libMesh;
  
  Biharmonic::JR::JR(EquationSystems& eqSys,
  		   const std::string& name_in,
  		   const unsigned int number_in) :
    TransientNonlinearImplicitSystem(eqSys,name_in,number_in),
    _biharmonic(dynamic_cast<Biharmonic&>(eqSys))
  {
  #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
    ERROR("Must have second derivatives enabled");
  #endif
  
  #ifdef LIBMESH_ENABLE_PERIODIC
    DofMap& dof_map = get_dof_map();
    PeriodicBoundary xbdry(RealVectorValue(1.0, 0.0, 0.0));
  #if LIBMESH_DIM > 1
    PeriodicBoundary ybdry(RealVectorValue(0.0, 1.0, 0.0));
  #endif
  #if LIBMESH_DIM > 2
    PeriodicBoundary zbdry(RealVectorValue(0.0, 0.0, 1.0));
  #endif
  
    switch(_biharmonic._dim)
      {
      case 1:
        xbdry.myboundary = 0;
        xbdry.pairedboundary = 1;
        dof_map.add_periodic_boundary(xbdry);
        break;
  #if LIBMESH_DIM > 1
      case 2:
        xbdry.myboundary = 3;
        xbdry.pairedboundary = 1;
        dof_map.add_periodic_boundary(xbdry);
        ybdry.myboundary = 0;
        ybdry.pairedboundary = 2;
        dof_map.add_periodic_boundary(ybdry);
        break;
  #endif
  #if LIBMESH_DIM > 2
      case 3:
        xbdry.myboundary = 4;
        xbdry.pairedboundary = 2;
        dof_map.add_periodic_boundary(xbdry);
        ybdry.myboundary = 1;
        ybdry.pairedboundary = 3;
        dof_map.add_periodic_boundary(ybdry);
        zbdry.myboundary = 0;
        zbdry.pairedboundary = 5;
        dof_map.add_periodic_boundary(zbdry);
        break;
  #endif
      default:
        libmesh_error();
      }
  #endif // LIBMESH_ENABLE_PERIODIC
  
  
    add_variable("u", THIRD, HERMITE);
  
    attach_init_object(*this);
  
    nonlinear_solver->residual_and_jacobian_object = this;
  
    nonlinear_solver->bounds_object = this;
  }
  
  
  
  
  
  void Biharmonic::JR::initialize()
  {
    if (_biharmonic._verbose)
      libMesh::out << ">>> Initializing Biharmonic::JR\n";
  
    Parameters parameters;
    parameters.set<Point>("center") = _biharmonic._initialCenter;
    parameters.set<Real>("width")   = _biharmonic._initialWidth;
  
    if (_biharmonic._initialState == Biharmonic::BALL)
      project_solution(Biharmonic::JR::InitialDensityBall, Biharmonic::JR::InitialGradientZero, parameters);
  
    if (_biharmonic._initialState == Biharmonic::ROD)
      project_solution(Biharmonic::JR::InitialDensityRod, Biharmonic::JR::InitialGradientZero, parameters);
  
    if (_biharmonic._initialState == Biharmonic::STRIP)
      project_solution(Biharmonic::JR::InitialDensityStrip, Biharmonic::JR::InitialGradientZero, parameters);
  
    *(old_local_solution) = *(current_local_solution);
  
    if (_biharmonic._verbose)
      libMesh::out << "<<< Initializing Biharmonic::JR\n";
  }
  
  
  
  
  
  
  Number Biharmonic::JR::InitialDensityBall(const Point& p,
  					  const Parameters& parameters,
  					  const std::string&,
  					  const std::string&)
  {
    Point center = parameters.get<Point>("center");
    Real width = parameters.get<Real>("width");
    Point pc = p-center;
    Real r = pc.size();
    return (r < width) ? 1.0 : -0.5;
  }
  
  
  
  
  Number Biharmonic::JR::InitialDensityRod(const Point& p,
  					 const Parameters& parameters,
  					 const std::string&,
  					 const std::string&)
  {
    Point center = parameters.get<Point>("center");
    Real width = parameters.get<Real>("width");
    Real r = sqrt((p(0)-center(0))*(p(0)-center(0)) + (p(1)-center(1))*(p(1)-center(1)));
    return (r < width) ? 1.0 : -0.5;
  }
  
  
  
  
  
  Number Biharmonic::JR::InitialDensityStrip(const Point& p,
  					   const Parameters& parameters,
  					   const std::string&,
  					   const std::string&)
  {
    Point center = parameters.get<Point>("center");
    Real width = parameters.get<Real>("width");
    Real r = sqrt((p(0)-center(0))*(p(0)-center(0)));
    return (r < width) ? 1.0 : -0.5;
  }
  
  
  
  
  Gradient Biharmonic::JR::InitialGradientZero(const Point&,
  					     const Parameters&,
  					     const std::string&,
  					     const std::string&)
  {
    return Gradient(0.0,0.0,0.0);
  }
  
  
  
  
  void Biharmonic::JR::residual_and_jacobian(const NumericVector<Number> &u,
  					   NumericVector<Number> *R,
  					   SparseMatrix<Number> *J,
  					   NonlinearImplicitSystem&)
  {
  #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
    if (!R && !J)
      return;
  
    PerfLog perf_log ("Biharmonic Residual and Jacobian", false);
  
    const DofMap& dof_map = get_dof_map();
  
    FEType fe_type = dof_map.variable_type(0);
  
    AutoPtr<FEBase> fe (FEBase::build(_biharmonic._dim, fe_type));
  
    AutoPtr<QBase> qrule(fe_type.default_quadrature_rule(_biharmonic._dim));
  
    fe->attach_quadrature_rule (qrule.get());
  
    const std::vector<Real>& JxW = fe->get_JxW();
  
    const std::vector<std::vector<Real> >& phi = fe->get_phi();
  
    const std::vector<std::vector<RealGradient> >& dphi = fe->get_dphi();
  
    const std::vector<std::vector<RealTensor> >& d2phi = fe->get_d2phi();
  
    std::vector<Real> Laplacian_phi_qp;
  
    DenseMatrix<Number> Je;
    DenseVector<Number> Re;
  
    std::vector<dof_id_type> dof_indices;
  
    const NumericVector<Number>& u_old = *old_local_solution;
  
  
    MeshBase::const_element_iterator       el     = _biharmonic._mesh->active_local_elements_begin();
    const MeshBase::const_element_iterator end_el = _biharmonic._mesh->active_local_elements_end();
  
    for ( ; el != end_el; ++el) {
      const Elem* elem = *el;
  
      dof_map.dof_indices (elem, dof_indices);
  
      fe->reinit (elem);
  
      if (J)
        Je.resize(dof_indices.size(), dof_indices.size());
  
      if (R)
        Re.resize(dof_indices.size());
  
      Laplacian_phi_qp.resize(dof_indices.size());
  
      for (unsigned int qp=0; qp<qrule->n_points(); qp++)
        {
  	Number u_qp = 0.0, u_old_qp = 0.0, Laplacian_u_qp = 0.0, Laplacian_u_old_qp = 0.0;
  	Gradient grad_u_qp(0.0,0.0,0.0), grad_u_old_qp(0.0,0.0,0.0);
  	Number M_qp = 1.0, M_old_qp = 1.0, M_prime_qp = 0.0, M_prime_old_qp = 0.0;
  
  	for (unsigned int i=0; i<phi.size(); i++)
  	  {
  	    Laplacian_phi_qp[i] = d2phi[i][qp](0,0);
  	    grad_u_qp(0) += u(dof_indices[i])*dphi[i][qp](0);
  	    grad_u_old_qp(0) += u_old(dof_indices[i])*dphi[i][qp](0);
  
  	    if (_biharmonic._dim > 1)
  	      {
  		Laplacian_phi_qp[i] += d2phi[i][qp](1,1);
  		grad_u_qp(1) += u(dof_indices[i])*dphi[i][qp](1);
  		grad_u_old_qp(1) += u_old(dof_indices[i])*dphi[i][qp](1);
  	      }
  	    if (_biharmonic._dim > 2)
  	      {
  		Laplacian_phi_qp[i] += d2phi[i][qp](2,2);
  		grad_u_qp(2) += u(dof_indices[i])*dphi[i][qp](2);
  		grad_u_old_qp(2) += u_old(dof_indices[i])*dphi[i][qp](2);
  	      }
  	    u_qp     += phi[i][qp]*u(dof_indices[i]);
  	    u_old_qp += phi[i][qp]*u_old(dof_indices[i]);
  	    Laplacian_u_qp     += Laplacian_phi_qp[i]*u(dof_indices[i]);
  	    Laplacian_u_old_qp += Laplacian_phi_qp[i]*u_old(dof_indices[i]);
  	  } // for i
  
  	if (_biharmonic._degenerate)
  	  {
  	    M_qp           = 1.0 - u_qp*u_qp;
  	    M_old_qp       = 1.0 - u_old_qp*u_old_qp;
  	    M_prime_qp     = -2.0*u_qp;
  	    M_prime_old_qp = -2.0*u_old_qp;
  	  }
  
  	for (unsigned int i=0; i<phi.size(); i++)
  	  {
  	    if (R)
  	      {
  		Number ri = 0.0, ri_old = 0.0;
  		ri     -= Laplacian_phi_qp[i]*M_qp*_biharmonic._kappa*Laplacian_u_qp;
  		ri_old -= Laplacian_phi_qp[i]*M_old_qp*_biharmonic._kappa*Laplacian_u_old_qp;
  
  		if (_biharmonic._degenerate)
  		  {
  		    ri       -= (dphi[i][qp]*grad_u_qp)*M_prime_qp*(_biharmonic._kappa*Laplacian_u_qp);
  		    ri_old   -= (dphi[i][qp]*grad_u_old_qp)*M_prime_old_qp*(_biharmonic._kappa*Laplacian_u_old_qp);
  		  }
  
  		if (_biharmonic._cahn_hillard)
  		  {
  		    if (_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
  		      {
  			ri += Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp;
  			ri_old += Laplacian_phi_qp[i]*M_old_qp*_biharmonic._theta_c*(u_old_qp*u_old_qp - 1.0)*u_old_qp;
  			if (_biharmonic._degenerate)
  			  {
  			    ri     += (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp;
  			    ri_old += (dphi[i][qp]*grad_u_old_qp)*M_prime_old_qp*_biharmonic._theta_c*(u_old_qp*u_old_qp - 1.0)*u_old_qp;
  			  }
  		      }// if(_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
  
  		    if (_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
  		      {
  			ri -= Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*u_qp;
  			ri_old -= Laplacian_phi_qp[i]*M_old_qp*_biharmonic._theta_c*u_old_qp;
  			if (_biharmonic._degenerate)
  			  {
  			    ri     -= (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*u_qp;
  			    ri_old -= (dphi[i][qp]*grad_u_old_qp)*M_prime_old_qp*_biharmonic._theta_c*u_old_qp;
  			  }
  		      } // if(_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
  
  		    if (_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
  		      {
  			switch(_biharmonic._log_truncation)
  			  {
  			  case 2:
  			    break;
  			  case 3:
  			    break;
  			  default:
  			    break;
  			  }// switch(_biharmonic._log_truncation)
  		      }// if(_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
  		  }// if(_biharmonic._cahn_hillard)
  		Re(i) += JxW[qp]*((u_qp-u_old_qp)*phi[i][qp]-_biharmonic._dt*0.5*((2.0-_biharmonic._cnWeight)*ri + _biharmonic._cnWeight*ri_old));
  	      } // if (R)
  
  	    if (J)
  	      {
  		Number M_prime_prime_qp = 0.0;
  		if(_biharmonic._degenerate) M_prime_prime_qp = -2.0;
  		for (unsigned int j=0; j<phi.size(); j++)
  		  {
  		    Number ri_j = 0.0;
  		    ri_j -= Laplacian_phi_qp[i]*M_qp*_biharmonic._kappa*Laplacian_phi_qp[j];
  		    if (_biharmonic._degenerate)
  		      {
  			ri_j -=
  			  Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._kappa*Laplacian_u_qp               +
  			  (dphi[i][qp]*dphi[j][qp])*M_prime_qp*(_biharmonic._kappa*Laplacian_u_qp)                  +
  			  (dphi[i][qp]*grad_u_qp)*(M_prime_prime_qp*phi[j][qp])*(_biharmonic._kappa*Laplacian_u_qp) +
  			  (dphi[i][qp]*grad_u_qp)*(M_prime_qp)*(_biharmonic._kappa*Laplacian_phi_qp[j]);
  		      }
  
  		    if (_biharmonic._cahn_hillard)
  		      {
  			if(_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
  			  {
  			    ri_j +=
  			      Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp +
  			      Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*(3.0*u_qp*u_qp - 1.0)*phi[j][qp]        +
  			      (dphi[i][qp]*dphi[j][qp])*M_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp      +
  			      (dphi[i][qp]*grad_u_qp)*M_prime_prime_qp*_biharmonic._theta_c*(u_qp*u_qp - 1.0)*u_qp  +
  			      (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*(3.0*u_qp*u_qp - 1.0)*phi[j][qp];
  			  }// if(_biharmonic._energy == DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_WELL)
  
  			if (_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
  			  {
  			    ri_j -=
  			      Laplacian_phi_qp[i]*M_prime_qp*phi[j][qp]*_biharmonic._theta_c*u_qp                   +
  			      Laplacian_phi_qp[i]*M_qp*_biharmonic._theta_c*phi[j][qp]                              +
  			      (dphi[i][qp]*dphi[j][qp])*M_prime_qp*_biharmonic._theta_c*u_qp                        +
  			      (dphi[i][qp]*grad_u_qp)*M_prime_prime_qp*_biharmonic._theta_c*u_qp                    +
  			      (dphi[i][qp]*grad_u_qp)*M_prime_qp*_biharmonic._theta_c*phi[j][qp];
  			  } // if(_biharmonic._energy == DOUBLE_OBSTACLE || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
  
  			if (_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
  			  {
  			    switch(_biharmonic._log_truncation)
  			      {
  			      case 2:
  				break;
  			      case 3:
  				break;
  			      default:
  				break;
  			      }// switch(_biharmonic._log_truncation)
  			  }// if(_biharmonic._energy == LOG_DOUBLE_WELL || _biharmonic._energy == LOG_DOUBLE_OBSTACLE)
  		      }// if(_biharmonic._cahn_hillard)
  		    Je(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp] - 0.5*_biharmonic._dt*(2.0-_biharmonic._cnWeight)*ri_j);
  		  } // for j
  	      } // if (J)
  	  } // for i
        } // for qp
  
      if (R)
        {
  	dof_map.constrain_element_vector(Re, dof_indices);
  	R->add_vector(Re, dof_indices);
        }
  
      if (J)
        {
  	dof_map.constrain_element_matrix(Je, dof_indices);
  	J->add_matrix(Je, dof_indices);
        }
    } // for el
  #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
  }
  
  
  
  
  
  void Biharmonic::JR::bounds(NumericVector<Number> &XL, NumericVector<Number>& XU, NonlinearImplicitSystem&)
  {
  
    PerfLog perf_log ("Biharmonic bounds", false);
  
    const DofMap& dof_map = get_dof_map();
  
    FEType fe_type = dof_map.variable_type(0);
  
    AutoPtr<FEBase> fe (FEBase::build(_biharmonic._dim, fe_type));
  
    DenseVector<Number> XLe, XUe;
  
    std::vector<dof_id_type> dof_indices;
  
    MeshBase::const_element_iterator       el     = _biharmonic._mesh->active_local_elements_begin();
    const MeshBase::const_element_iterator end_el = _biharmonic._mesh->active_local_elements_end();
  
    for ( ; el != end_el; ++el)
      {
        const std::vector<std::vector<Real> >& phi = fe->get_phi();
  
        dof_map.dof_indices (*el, dof_indices);
  
        XLe.resize(dof_indices.size());
        XUe.resize(dof_indices.size());
  
        std::vector<Point> nodes;
        fe->get_refspace_nodes((*el)->type(), nodes);
  
        fe->reinit(*el, &nodes);
  
        for (unsigned int i=0; i<phi.size(); ++i)
  	{
  	  Real infinity = 1.0e20;
  	  Real bound = infinity;
  	  for(unsigned int j = 0; j < nodes.size(); ++j) {
  	    if(phi[i][j]) {
  	      bound = 1.0/fabs(phi[i][j]);
  	      break;
  	    }
  	  }
  
  	  XLe(i) = -bound;
  	  XUe(i) = bound;
  	}
        XL.insert(XLe, dof_indices);
        XU.insert(XUe, dof_indices);
      }
  }



The source file miscellaneous_ex7.C without comments:

 
  #include "biharmonic.h"
  
  
  using namespace libMesh;
  
  void print_help(int argc, char** argv);
  
  int main(int argc, char** argv)
  {
    LibMeshInit init (argc, argv);
    if (on_command_line("--help"))
      print_help(argc, argv);
    else
      {
  #if !defined(LIBMESH_ENABLE_SECOND_DERIVATIVES)
        libmesh_example_assert(false, "--enable-second");
  #elif !defined(LIBMESH_ENABLE_PERIODIC)
        libmesh_example_assert(false, "--enable-periodic");
  #endif
  
        libmesh_example_assert(libMesh::default_solver_package() == PETSC_SOLVERS, "--enable-petsc");
  
        const int dim = command_line_value("dim",1);
  
        libmesh_example_assert(dim <= LIBMESH_DIM, "2D/3D support");
  
        Biharmonic* biharmonic;
        Biharmonic::Create(&biharmonic,init.comm());
        biharmonic->viewParameters();
        biharmonic->init();
        biharmonic->run();
        Biharmonic::Destroy(&biharmonic);
      }
    return 0;
  }
  
  
  
  
  
  void print_help(int, char** argv)
  {
    libMesh::out << "This example solves the Cahn-Hillard equation with chemical potential f:\n"
  	       << "    u_t = \\div(M(u)\\grad f(u))\n"
  	       << "Here we have\n"
  	       << "    u, -1 <= u <= 1        -- relative concentration (difference of two concentrations in a binary mixture) \n"
  	       << "    M, M >= 0              -- mobility of the mixture\n"
  	       << "    f = \\delta E/\\delta u  -- variational derivative of the free energy functional E\n"
  	       << "    E = \\int[\\kappa/2 |\\grac u|^ + g(u)]\n"
  	       << "where the gradient term is the interfacial energy density with \\kappa quantifying the energy of the interface,\n"
  	       << "and g(u) is the bulk energy density\n"
  	       << "    g(u) = \\theta L(u) + \\theta_c W(u),\n"
  	       << "L(u) is the (optional, in this model) logarithmic term corresponding to the entropy of the mixture:\n"
  	       << "    L(u) = (\\theta/2)[(1+u)\\ln((1+u)/2) + (1-u)\\ln((1-u)/2)],\n"
  	       << "where \\theta is related to the Boltzmann factor k_B T - a proxy for the absolute temperature T.\n"
  	       << "L can be optionally approximated ('truncated') using a quadratic or a cubic polynomial on [-1,1]\n"
  	       << "W(u) is the (optional, in this model) potential promoting demixing.  It can take the form of \n"
  	       << "a 'double well' potential\n"
  	       << "    W(u) = \\theta_c (u^4/4 - u^2/2),\n"
  	       << "         or \n"
  	       << "a 'double obstacle' potential\n"
  	       << "    W(u) = (\\theta_c/2)(1-u^2),\n"
  	       << "where \\theta_c is the critical 'temperature'.\n"
  	       << "Finally, mobility M can be constant of 'degenerate', by which we mean that M is varying with u and \n"
  	       << "vanishing (degenerating) whenever u reaches critical values +1 or -1:\n"
  	       << "    M(u) = 1.0\n"
  	       << "      or\n"
  	       << "    M(u) = (1.0 - u^2)\n"
  	       << "Degenerate mobility should generally be used only in conjunction with logarithmic free energy terms.\n\n"
  	       << "The equation is solved on a periodic domain (in 1D, 2D or 3D)\n"
  	       << "using a Galerkin formulation with C^1 elements approximating the H^2_{per} function space.\n\n"
  	       << "\n-----------\n"
  	       << "COMPILING: "
  	       << "\n-----------\n"
  	       << "Compile as follows (assuming libmesh has been built): \n"
  	       << "METHOD=<method> make \n"
  	       << "where <method> is dbg or opt.\n"
  	       << "\n-----------\n"
  	       << "HELP:        "
  	       << "\n-----------\n"
  	       << "Print this help message:\n"
  	       << argv[0] << " --help\n"
  	       << "\n-----------\n"
  	       << "RUNNING:     "
  	       << "\n-----------\n"
  	       << "Run in serial with build METHOD <method> as follows:\n"
  	       << "\n"
  	       << argv[0] << "\n"
  	       << "               [--verbose] dim=<1|2|3> N=<number_of_linear_elements> \n"
  	       << "               kappa=<kappa_value> growth=<yes|no> degenerate=<yes|no> [--cahn-hillard]                                           \n"
  	       << "               [--netforce]  energy=<double_well|double_obstacle|log_double_well|log_double_obstacle>  log_truncation_order=<2|3> \n"
  	       << "               theta=<theta_value> theta_c=<theta_c_value>                                                                        \n"
  	       << "               initial_state=<ball|rod|strip> initial_center='x [y [z]]' initial_width=<width>                                    \n"
  	       << "               min_time=<initial_time> max_time=<final_time> dt=<timestep_size> crank_nicholson_weight=<between_0_and_1>          \n"
  	       << "               output_base=<base_filename> output_dt=<output_timestep_size> [--use-petsc-dm -snes_type virs]                      \n"
  	       << "\n"
  	       << argv[0] << " --verbose \n"
  	       << "is a pretty good start.\n"
  	       << "\nModeling a 1D system with 2D or 3D (for a strip the second and third components of the center are immaterial):\n"
  	       << argv[0]<< " --verbose dim=1 N=1024 initial_state=strip initial_center=0.5 initial_width=0.1 dt=1e-10 max_time=1e-6\n"
  	       << argv[0]<< " --verbose dim=2 N=64   initial_state=strip initial_center=0.5 initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
  	       << argv[0]<< " --verbose dim=3 N=32   initial_state=strip initial_center=0.5 initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
  	       << "\n"
  	       << "Modeling a 2D system with 3D (for a rod the third component of the center is immaterial) \n"
  	       << argv[0]<< " --verbose dim=2 N=64   initial_state=rod initial_center='0.5 0.5' initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
  	       << argv[0]<< " --verbose dim=3 N=32   initial_state=rod initial_center='0.5 0.5' initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
  	       << "\n"
  	       << "A 3D system with an initial ball in the center\n"
  	       << argv[0] << " --verbose dim=3 N=32   initial_state=ball initial_center='0.5 0.5 0.5' initial_width=0.1 dt=1e-10 max_time=1e-6 \n"
  	       << "\n"
  	       << "Add --use-petsc-dm -snes_type virs to run the variational inequality version that ensures the solution is between -1.0 and 1.0 at all times.\n\n"
  	       << std::endl;
  }



The console output of the program:

make[4]: Entering directory `/net/spark/workspace/roystgnr/libmesh/git/devel/examples/miscellaneous/miscellaneous_ex7'
***************************************************************
* Running Example miscellaneous_ex7:
*  mpirun -np 4 example-devel --verbose dim=1 N=1024 initial_state=strip initial_center=0.5 initial_width=0.1 dt=1e-10 max_time=1e-8 -pc_type bjacobi -sub_pc_type ilu -sub_pc_factor_levels 4 -sub_pc_factor_zeropivot 0 -ksp_right_pc
***************************************************************
 
Biharmonic parameters:
verbose mode is on
mesh dimension           = 1
initial linear mesh size = 1024
kappa                    = 1
growth                   = 0
degenerate               = 0
Cahn-Hillard             = 0
netforce                 = 0
energy                   = 1
tol                      = 1e-08
theta                    = 0.001
theta_c                  = 1
log truncation           = 2
initial timestep size    = 1e-10
initial state:             strip
initial state center     = 0.5
initial state width      = 0.1
initial time (min_time)  = 0
integration time         = 1e-08
final time   (max_time)  = 1e-08
Crank-Nicholson weight   = 1
Output timestep          = 0
Output filename base:      bih.1
>>> Initializing Biharmonic
>>> Initializing Biharmonic::JR
<<< Initializing Biharmonic::JR
<<< Initializing Biharmonic
Writing state 1 at time 0 to file bih.1.e; output a total of 1 states so far
Solving for state 2, time 0
  NL step  0, |residual|_2 = 3.865471e+00
  NL step  1, |residual|_2 = 6.791353e-15
Writing state 2 at time 1.000000e-10 to file bih.1.e; output a total of 2 states so far
Solving for state 3, time 1.000000e-10
  NL step  0, |residual|_2 = 3.861870e+00
  NL step  1, |residual|_2 = 5.295328e-13
Writing state 3 at time 2.000000e-10 to file bih.1.e; output a total of 3 states so far
Solving for state 4, time 2.000000e-10
  NL step  0, |residual|_2 = 3.858634e+00
  NL step  1, |residual|_2 = 6.826469e-12
Writing state 4 at time 3.000000e-10 to file bih.1.e; output a total of 4 states so far
Solving for state 5, time 3.000000e-10
  NL step  0, |residual|_2 = 3.855483e+00
  NL step  1, |residual|_2 = 2.625789e-11
Writing state 5 at time 4.000000e-10 to file bih.1.e; output a total of 5 states so far
Solving for state 6, time 4.000000e-10
  NL step  0, |residual|_2 = 3.852464e+00
  NL step  1, |residual|_2 = 3.298126e-11
Writing state 6 at time 5.000000e-10 to file bih.1.e; output a total of 6 states so far
Solving for state 7, time 5.000000e-10
  NL step  0, |residual|_2 = 3.849501e+00
  NL step  1, |residual|_2 = 5.351723e-10
Writing state 7 at time 6.000000e-10 to file bih.1.e; output a total of 7 states so far
Solving for state 8, time 6.000000e-10
  NL step  0, |residual|_2 = 3.846621e+00
  NL step  1, |residual|_2 = 1.678735e-09
Writing state 8 at time 7.000000e-10 to file bih.1.e; output a total of 8 states so far
Solving for state 9, time 7.000000e-10
  NL step  0, |residual|_2 = 3.843784e+00
  NL step  1, |residual|_2 = 1.776743e-09
Writing state 9 at time 8.000000e-10 to file bih.1.e; output a total of 9 states so far
Solving for state 10, time 8.000000e-10
  NL step  0, |residual|_2 = 3.841009e+00
  NL step  1, |residual|_2 = 2.711450e-09
Writing state 10 at time 9.000000e-10 to file bih.1.e; output a total of 10 states so far
Solving for state 11, time 9.000000e-10
  NL step  0, |residual|_2 = 3.838268e+00
  NL step  1, |residual|_2 = 1.204330e-08
Writing state 11 at time 1.000000e-09 to file bih.1.e; output a total of 11 states so far
Solving for state 12, time 1.000000e-09
  NL step  0, |residual|_2 = 3.835577e+00
  NL step  1, |residual|_2 = 1.982070e-08
Writing state 12 at time 1.100000e-09 to file bih.1.e; output a total of 12 states so far
Solving for state 13, time 1.100000e-09
  NL step  0, |residual|_2 = 3.832915e+00
  NL step  1, |residual|_2 = 1.887360e-08
Writing state 13 at time 1.200000e-09 to file bih.1.e; output a total of 13 states so far
Solving for state 14, time 1.200000e-09
  NL step  0, |residual|_2 = 3.830295e+00
  NL step  1, |residual|_2 = 7.909848e-09
Writing state 14 at time 1.300000e-09 to file bih.1.e; output a total of 14 states so far
Solving for state 15, time 1.300000e-09
  NL step  0, |residual|_2 = 3.827700e+00
  NL step  1, |residual|_2 = 1.049956e-08
Writing state 15 at time 1.400000e-09 to file bih.1.e; output a total of 15 states so far
Solving for state 16, time 1.400000e-09
  NL step  0, |residual|_2 = 3.825140e+00
  NL step  1, |residual|_2 = 3.280321e-08
Writing state 16 at time 1.500000e-09 to file bih.1.e; output a total of 16 states so far
Solving for state 17, time 1.500000e-09
  NL step  0, |residual|_2 = 3.822604e+00
  NL step  1, |residual|_2 = 5.475702e-08
  NL step  2, |residual|_2 = 2.323455e-15
Writing state 17 at time 1.600000e-09 to file bih.1.e; output a total of 17 states so far
Solving for state 18, time 1.600000e-09
  NL step  0, |residual|_2 = 3.820098e+00
  NL step  1, |residual|_2 = 3.707085e-08
Writing state 18 at time 1.700000e-09 to file bih.1.e; output a total of 18 states so far
Solving for state 19, time 1.700000e-09
  NL step  0, |residual|_2 = 3.817613e+00
  NL step  1, |residual|_2 = 9.270300e-08
  NL step  2, |residual|_2 = 2.353413e-15
Writing state 19 at time 1.800000e-09 to file bih.1.e; output a total of 19 states so far
Solving for state 20, time 1.800000e-09
  NL step  0, |residual|_2 = 3.815156e+00
  NL step  1, |residual|_2 = 1.778004e-08
Writing state 20 at time 1.900000e-09 to file bih.1.e; output a total of 20 states so far
Solving for state 21, time 1.900000e-09
  NL step  0, |residual|_2 = 3.812718e+00
  NL step  1, |residual|_2 = 1.166306e-07
  NL step  2, |residual|_2 = 2.344579e-15
Writing state 21 at time 2.000000e-09 to file bih.1.e; output a total of 21 states so far
Solving for state 22, time 2.000000e-09
  NL step  0, |residual|_2 = 3.810305e+00
  NL step  1, |residual|_2 = 5.194027e-08
  NL step  2, |residual|_2 = 2.292010e-15
Writing state 22 at time 2.100000e-09 to file bih.1.e; output a total of 22 states so far
Solving for state 23, time 2.100000e-09
  NL step  0, |residual|_2 = 3.807909e+00
  NL step  1, |residual|_2 = 7.569577e-08
  NL step  2, |residual|_2 = 2.495208e-15
Writing state 23 at time 2.200000e-09 to file bih.1.e; output a total of 23 states so far
Solving for state 24, time 2.200000e-09
  NL step  0, |residual|_2 = 3.805536e+00
  NL step  1, |residual|_2 = 7.326139e-08
  NL step  2, |residual|_2 = 2.341764e-15
Writing state 24 at time 2.300000e-09 to file bih.1.e; output a total of 24 states so far
Solving for state 25, time 2.300000e-09
  NL step  0, |residual|_2 = 3.803179e+00
  NL step  1, |residual|_2 = 4.221623e-08
  NL step  2, |residual|_2 = 2.358325e-15
Writing state 25 at time 2.400000e-09 to file bih.1.e; output a total of 25 states so far
Solving for state 26, time 2.400000e-09
  NL step  0, |residual|_2 = 3.800842e+00
  NL step  1, |residual|_2 = 1.024370e-07
  NL step  2, |residual|_2 = 2.288452e-15
Writing state 26 at time 2.500000e-09 to file bih.1.e; output a total of 26 states so far
Solving for state 27, time 2.500000e-09
  NL step  0, |residual|_2 = 3.798522e+00
  NL step  1, |residual|_2 = 3.410027e-08
Writing state 27 at time 2.600000e-09 to file bih.1.e; output a total of 27 states so far
Solving for state 28, time 2.600000e-09
  NL step  0, |residual|_2 = 3.796220e+00
  NL step  1, |residual|_2 = 2.482465e-07
  NL step  2, |residual|_2 = 2.646697e-15
Writing state 28 at time 2.700000e-09 to file bih.1.e; output a total of 28 states so far
Solving for state 29, time 2.700000e-09
  NL step  0, |residual|_2 = 3.793932e+00
  NL step  1, |residual|_2 = 1.025109e-07
  NL step  2, |residual|_2 = 2.243753e-15
Writing state 29 at time 2.800000e-09 to file bih.1.e; output a total of 29 states so far
Solving for state 30, time 2.800000e-09
  NL step  0, |residual|_2 = 3.791663e+00
  NL step  1, |residual|_2 = 2.505431e-07
  NL step  2, |residual|_2 = 2.438573e-15
Writing state 30 at time 2.900000e-09 to file bih.1.e; output a total of 30 states so far
Solving for state 31, time 2.900000e-09
  NL step  0, |residual|_2 = 3.789407e+00
  NL step  1, |residual|_2 = 9.550927e-08
  NL step  2, |residual|_2 = 2.429171e-15
Writing state 31 at time 3.000000e-09 to file bih.1.e; output a total of 31 states so far
Solving for state 32, time 3.000000e-09
  NL step  0, |residual|_2 = 3.787167e+00
  NL step  1, |residual|_2 = 2.581924e-07
  NL step  2, |residual|_2 = 2.508660e-15
Writing state 32 at time 3.100000e-09 to file bih.1.e; output a total of 32 states so far
Solving for state 33, time 3.100000e-09
  NL step  0, |residual|_2 = 3.784941e+00
  NL step  1, |residual|_2 = 9.692862e-08
  NL step  2, |residual|_2 = 2.404583e-15
Writing state 33 at time 3.200000e-09 to file bih.1.e; output a total of 33 states so far
Solving for state 34, time 3.200000e-09
  NL step  0, |residual|_2 = 3.782730e+00
  NL step  1, |residual|_2 = 2.594179e-07
  NL step  2, |residual|_2 = 2.491096e-15
Writing state 34 at time 3.300000e-09 to file bih.1.e; output a total of 34 states so far
Solving for state 35, time 3.300000e-09
  NL step  0, |residual|_2 = 3.780531e+00
  NL step  1, |residual|_2 = 9.887434e-08
  NL step  2, |residual|_2 = 2.479266e-15
Writing state 35 at time 3.400000e-09 to file bih.1.e; output a total of 35 states so far
Solving for state 36, time 3.400000e-09
  NL step  0, |residual|_2 = 3.778347e+00
  NL step  1, |residual|_2 = 2.527746e-07
  NL step  2, |residual|_2 = 2.569841e-15
Writing state 36 at time 3.500000e-09 to file bih.1.e; output a total of 36 states so far
Solving for state 37, time 3.500000e-09
  NL step  0, |residual|_2 = 3.776174e+00
  NL step  1, |residual|_2 = 9.893108e-08
  NL step  2, |residual|_2 = 2.460867e-15
Writing state 37 at time 3.600000e-09 to file bih.1.e; output a total of 37 states so far
Solving for state 38, time 3.600000e-09
  NL step  0, |residual|_2 = 3.774015e+00
  NL step  1, |residual|_2 = 2.386953e-07
  NL step  2, |residual|_2 = 2.287197e-15
Writing state 38 at time 3.700000e-09 to file bih.1.e; output a total of 38 states so far
Solving for state 39, time 3.700000e-09
  NL step  0, |residual|_2 = 3.771868e+00
  NL step  1, |residual|_2 = 9.823210e-08
  NL step  2, |residual|_2 = 2.425269e-15
Writing state 39 at time 3.800000e-09 to file bih.1.e; output a total of 39 states so far
Solving for state 40, time 3.800000e-09
  NL step  0, |residual|_2 = 3.769733e+00
  NL step  1, |residual|_2 = 2.184248e-07
  NL step  2, |residual|_2 = 2.412187e-15
Writing state 40 at time 3.900000e-09 to file bih.1.e; output a total of 40 states so far
Solving for state 41, time 3.900000e-09
  NL step  0, |residual|_2 = 3.767609e+00
  NL step  1, |residual|_2 = 9.982364e-08
  NL step  2, |residual|_2 = 2.594741e-15
Writing state 41 at time 4.000000e-09 to file bih.1.e; output a total of 41 states so far
Solving for state 42, time 4.000000e-09
  NL step  0, |residual|_2 = 3.765498e+00
  NL step  1, |residual|_2 = 1.936364e-07
  NL step  2, |residual|_2 = 2.566911e-15
Writing state 42 at time 4.100000e-09 to file bih.1.e; output a total of 42 states so far
Solving for state 43, time 4.100000e-09
  NL step  0, |residual|_2 = 3.763397e+00
  NL step  1, |residual|_2 = 1.069263e-07
  NL step  2, |residual|_2 = 2.410036e-15
Writing state 43 at time 4.200000e-09 to file bih.1.e; output a total of 43 states so far
Solving for state 44, time 4.200000e-09
  NL step  0, |residual|_2 = 3.761307e+00
  NL step  1, |residual|_2 = 1.662948e-07
  NL step  2, |residual|_2 = 2.717666e-15
Writing state 44 at time 4.300000e-09 to file bih.1.e; output a total of 44 states so far
Solving for state 45, time 4.300000e-09
  NL step  0, |residual|_2 = 3.759228e+00
  NL step  1, |residual|_2 = 1.210237e-07
  NL step  2, |residual|_2 = 2.503640e-15
Writing state 45 at time 4.400000e-09 to file bih.1.e; output a total of 45 states so far
Solving for state 46, time 4.400000e-09
  NL step  0, |residual|_2 = 3.757160e+00
  NL step  1, |residual|_2 = 1.387212e-07
  NL step  2, |residual|_2 = 2.443778e-15
Writing state 46 at time 4.500000e-09 to file bih.1.e; output a total of 46 states so far
Solving for state 47, time 4.500000e-09
  NL step  0, |residual|_2 = 3.755101e+00
  NL step  1, |residual|_2 = 1.413200e-07
  NL step  2, |residual|_2 = 2.453897e-15
Writing state 47 at time 4.600000e-09 to file bih.1.e; output a total of 47 states so far
Solving for state 48, time 4.600000e-09
  NL step  0, |residual|_2 = 3.753053e+00
  NL step  1, |residual|_2 = 1.138969e-07
  NL step  2, |residual|_2 = 2.520689e-15
Writing state 48 at time 4.700000e-09 to file bih.1.e; output a total of 48 states so far
Solving for state 49, time 4.700000e-09
  NL step  0, |residual|_2 = 3.751014e+00
  NL step  1, |residual|_2 = 1.658550e-07
  NL step  2, |residual|_2 = 2.425168e-15
Writing state 49 at time 4.800000e-09 to file bih.1.e; output a total of 49 states so far
Solving for state 50, time 4.800000e-09
  NL step  0, |residual|_2 = 3.748985e+00
  NL step  1, |residual|_2 = 9.588393e-08
  NL step  2, |residual|_2 = 2.380361e-15
Writing state 50 at time 4.900000e-09 to file bih.1.e; output a total of 50 states so far
Solving for state 51, time 4.900000e-09
  NL step  0, |residual|_2 = 3.746966e+00
  NL step  1, |residual|_2 = 1.926354e-07
  NL step  2, |residual|_2 = 2.380779e-15
Writing state 51 at time 5.000000e-09 to file bih.1.e; output a total of 51 states so far
Solving for state 52, time 5.000000e-09
  NL step  0, |residual|_2 = 3.744956e+00
  NL step  1, |residual|_2 = 8.914698e-08
  NL step  2, |residual|_2 = 2.400894e-15
Writing state 52 at time 5.100000e-09 to file bih.1.e; output a total of 52 states so far
Solving for state 53, time 5.100000e-09
  NL step  0, |residual|_2 = 3.742955e+00
  NL step  1, |residual|_2 = 2.200581e-07
  NL step  2, |residual|_2 = 2.622921e-15
Writing state 53 at time 5.200000e-09 to file bih.1.e; output a total of 53 states so far
Solving for state 54, time 5.200000e-09
  NL step  0, |residual|_2 = 3.740963e+00
  NL step  1, |residual|_2 = 9.501731e-08
  NL step  2, |residual|_2 = 2.551815e-15
Writing state 54 at time 5.300000e-09 to file bih.1.e; output a total of 54 states so far
Solving for state 55, time 5.300000e-09
  NL step  0, |residual|_2 = 3.738979e+00
  NL step  1, |residual|_2 = 2.469333e-07
  NL step  2, |residual|_2 = 2.529382e-15
Writing state 55 at time 5.400000e-09 to file bih.1.e; output a total of 55 states so far
Solving for state 56, time 5.400000e-09
  NL step  0, |residual|_2 = 3.737005e+00
  NL step  1, |residual|_2 = 1.098357e-07
  NL step  2, |residual|_2 = 2.611139e-15
Writing state 56 at time 5.500000e-09 to file bih.1.e; output a total of 56 states so far
Solving for state 57, time 5.500000e-09
  NL step  0, |residual|_2 = 3.735039e+00
  NL step  1, |residual|_2 = 2.724056e-07
  NL step  2, |residual|_2 = 2.504157e-15
Writing state 57 at time 5.600000e-09 to file bih.1.e; output a total of 57 states so far
Solving for state 58, time 5.600000e-09
  NL step  0, |residual|_2 = 3.733082e+00
  NL step  1, |residual|_2 = 1.287789e-07
  NL step  2, |residual|_2 = 2.315606e-15
Writing state 58 at time 5.700000e-09 to file bih.1.e; output a total of 58 states so far
Solving for state 59, time 5.700000e-09
  NL step  0, |residual|_2 = 3.731132e+00
  NL step  1, |residual|_2 = 2.958737e-07
  NL step  2, |residual|_2 = 2.431640e-15
Writing state 59 at time 5.800000e-09 to file bih.1.e; output a total of 59 states so far
Solving for state 60, time 5.800000e-09
  NL step  0, |residual|_2 = 3.729191e+00
  NL step  1, |residual|_2 = 1.486029e-07
  NL step  2, |residual|_2 = 2.465803e-15
Writing state 60 at time 5.900000e-09 to file bih.1.e; output a total of 60 states so far
Solving for state 61, time 5.900000e-09
  NL step  0, |residual|_2 = 3.727257e+00
  NL step  1, |residual|_2 = 3.169299e-07
  NL step  2, |residual|_2 = 2.622898e-15
Writing state 61 at time 6.000000e-09 to file bih.1.e; output a total of 61 states so far
Solving for state 62, time 6.000000e-09
  NL step  0, |residual|_2 = 3.725332e+00
  NL step  1, |residual|_2 = 1.675296e-07
  NL step  2, |residual|_2 = 2.414417e-15
Writing state 62 at time 6.100000e-09 to file bih.1.e; output a total of 62 states so far
Solving for state 63, time 6.100000e-09
  NL step  0, |residual|_2 = 3.723414e+00
  NL step  1, |residual|_2 = 3.353144e-07
  NL step  2, |residual|_2 = 2.456989e-15
Writing state 63 at time 6.200000e-09 to file bih.1.e; output a total of 63 states so far
Solving for state 64, time 6.200000e-09
  NL step  0, |residual|_2 = 3.721504e+00
  NL step  1, |residual|_2 = 1.846217e-07
  NL step  2, |residual|_2 = 2.613437e-15
Writing state 64 at time 6.300000e-09 to file bih.1.e; output a total of 64 states so far
Solving for state 65, time 6.300000e-09
  NL step  0, |residual|_2 = 3.719601e+00
  NL step  1, |residual|_2 = 3.508814e-07
  NL step  2, |residual|_2 = 2.482544e-15
Writing state 65 at time 6.400000e-09 to file bih.1.e; output a total of 65 states so far
Solving for state 66, time 6.400000e-09
  NL step  0, |residual|_2 = 3.717706e+00
  NL step  1, |residual|_2 = 1.993937e-07
  NL step  2, |residual|_2 = 2.434944e-15
Writing state 66 at time 6.500000e-09 to file bih.1.e; output a total of 66 states so far
Solving for state 67, time 6.500000e-09
  NL step  0, |residual|_2 = 3.715818e+00
  NL step  1, |residual|_2 = 3.635735e-07
  NL step  2, |residual|_2 = 2.556709e-15
Writing state 67 at time 6.600000e-09 to file bih.1.e; output a total of 67 states so far
Solving for state 68, time 6.600000e-09
  NL step  0, |residual|_2 = 3.713937e+00
  NL step  1, |residual|_2 = 2.116130e-07
  NL step  2, |residual|_2 = 2.473761e-15
Writing state 68 at time 6.700000e-09 to file bih.1.e; output a total of 68 states so far
Solving for state 69, time 6.700000e-09
  NL step  0, |residual|_2 = 3.712063e+00
  NL step  1, |residual|_2 = 3.734010e-07
  NL step  2, |residual|_2 = 2.527880e-15
Writing state 69 at time 6.800000e-09 to file bih.1.e; output a total of 69 states so far
Solving for state 70, time 6.800000e-09
  NL step  0, |residual|_2 = 3.710197e+00
  NL step  1, |residual|_2 = 2.211980e-07
  NL step  2, |residual|_2 = 2.400609e-15
Writing state 70 at time 6.900000e-09 to file bih.1.e; output a total of 70 states so far
Solving for state 71, time 6.900000e-09
  NL step  0, |residual|_2 = 3.708336e+00
  NL step  1, |residual|_2 = 3.804259e-07
  NL step  2, |residual|_2 = 2.394281e-15
Writing state 71 at time 7.000000e-09 to file bih.1.e; output a total of 71 states so far
Solving for state 72, time 7.000000e-09
  NL step  0, |residual|_2 = 3.706484e+00
  NL step  1, |residual|_2 = 2.281627e-07
  NL step  2, |residual|_2 = 2.563391e-15
Writing state 72 at time 7.100000e-09 to file bih.1.e; output a total of 72 states so far
Solving for state 73, time 7.100000e-09
  NL step  0, |residual|_2 = 3.704637e+00
  NL step  1, |residual|_2 = 3.847491e-07
  NL step  2, |residual|_2 = 2.486401e-15
Writing state 73 at time 7.200000e-09 to file bih.1.e; output a total of 73 states so far
Solving for state 74, time 7.200000e-09
  NL step  0, |residual|_2 = 3.702797e+00
  NL step  1, |residual|_2 = 2.325837e-07
  NL step  2, |residual|_2 = 2.406717e-15
Writing state 74 at time 7.300000e-09 to file bih.1.e; output a total of 74 states so far
Solving for state 75, time 7.300000e-09
  NL step  0, |residual|_2 = 3.700963e+00
  NL step  1, |residual|_2 = 3.864999e-07
  NL step  2, |residual|_2 = 2.434656e-15
Writing state 75 at time 7.400000e-09 to file bih.1.e; output a total of 75 states so far
Solving for state 76, time 7.400000e-09
  NL step  0, |residual|_2 = 3.699137e+00
  NL step  1, |residual|_2 = 2.345800e-07
  NL step  2, |residual|_2 = 2.600773e-15
Writing state 76 at time 7.500000e-09 to file bih.1.e; output a total of 76 states so far
Solving for state 77, time 7.500000e-09
  NL step  0, |residual|_2 = 3.697316e+00
  NL step  1, |residual|_2 = 3.858276e-07
  NL step  2, |residual|_2 = 2.479296e-15
Writing state 77 at time 7.600000e-09 to file bih.1.e; output a total of 77 states so far
Solving for state 78, time 7.600000e-09
  NL step  0, |residual|_2 = 3.695502e+00
  NL step  1, |residual|_2 = 2.342992e-07
  NL step  2, |residual|_2 = 2.388925e-15
Writing state 78 at time 7.700000e-09 to file bih.1.e; output a total of 78 states so far
Solving for state 79, time 7.700000e-09
  NL step  0, |residual|_2 = 3.693694e+00
  NL step  1, |residual|_2 = 3.828952e-07
  NL step  2, |residual|_2 = 2.600408e-15
Writing state 79 at time 7.800000e-09 to file bih.1.e; output a total of 79 states so far
Solving for state 80, time 7.800000e-09
  NL step  0, |residual|_2 = 3.691892e+00
  NL step  1, |residual|_2 = 2.319080e-07
  NL step  2, |residual|_2 = 2.659483e-15
Writing state 80 at time 7.900000e-09 to file bih.1.e; output a total of 80 states so far
Solving for state 81, time 7.900000e-09
  NL step  0, |residual|_2 = 3.690096e+00
  NL step  1, |residual|_2 = 3.778738e-07
  NL step  2, |residual|_2 = 2.615159e-15
Writing state 81 at time 8.000000e-09 to file bih.1.e; output a total of 81 states so far
Solving for state 82, time 8.000000e-09
  NL step  0, |residual|_2 = 3.688306e+00
  NL step  1, |residual|_2 = 2.275868e-07
  NL step  2, |residual|_2 = 2.439835e-15
Writing state 82 at time 8.100000e-09 to file bih.1.e; output a total of 82 states so far
Solving for state 83, time 8.100000e-09
  NL step  0, |residual|_2 = 3.686522e+00
  NL step  1, |residual|_2 = 3.709386e-07
  NL step  2, |residual|_2 = 2.363289e-15
Writing state 83 at time 8.200000e-09 to file bih.1.e; output a total of 83 states so far
Solving for state 84, time 8.200000e-09
  NL step  0, |residual|_2 = 3.684744e+00
  NL step  1, |residual|_2 = 2.215246e-07
  NL step  2, |residual|_2 = 2.408869e-15
Writing state 84 at time 8.300000e-09 to file bih.1.e; output a total of 84 states so far
Solving for state 85, time 8.300000e-09
  NL step  0, |residual|_2 = 3.682971e+00
  NL step  1, |residual|_2 = 3.622657e-07
  NL step  2, |residual|_2 = 2.323864e-15
Writing state 85 at time 8.400000e-09 to file bih.1.e; output a total of 85 states so far
Solving for state 86, time 8.400000e-09
  NL step  0, |residual|_2 = 3.681205e+00
  NL step  1, |residual|_2 = 2.139174e-07
  NL step  2, |residual|_2 = 2.295697e-15
Writing state 86 at time 8.500000e-09 to file bih.1.e; output a total of 86 states so far
Solving for state 87, time 8.500000e-09
  NL step  0, |residual|_2 = 3.679443e+00
  NL step  1, |residual|_2 = 3.520296e-07
  NL step  2, |residual|_2 = 2.376511e-15
Writing state 87 at time 8.600000e-09 to file bih.1.e; output a total of 87 states so far
Solving for state 88, time 8.600000e-09
  NL step  0, |residual|_2 = 3.677688e+00
  NL step  1, |residual|_2 = 2.049672e-07
  NL step  2, |residual|_2 = 2.508000e-15
Writing state 88 at time 8.700000e-09 to file bih.1.e; output a total of 88 states so far
Solving for state 89, time 8.700000e-09
  NL step  0, |residual|_2 = 3.675938e+00
  NL step  1, |residual|_2 = 3.404019e-07
  NL step  2, |residual|_2 = 2.459481e-15
Writing state 89 at time 8.800000e-09 to file bih.1.e; output a total of 89 states so far
Solving for state 90, time 8.800000e-09
  NL step  0, |residual|_2 = 3.674193e+00
  NL step  1, |residual|_2 = 1.948828e-07
  NL step  2, |residual|_2 = 2.472973e-15
Writing state 90 at time 8.900000e-09 to file bih.1.e; output a total of 90 states so far
Solving for state 91, time 8.900000e-09
  NL step  0, |residual|_2 = 3.672454e+00
  NL step  1, |residual|_2 = 3.275496e-07
  NL step  2, |residual|_2 = 2.655535e-15
Writing state 91 at time 9.000000e-09 to file bih.1.e; output a total of 91 states so far
Solving for state 92, time 9.000000e-09
  NL step  0, |residual|_2 = 3.670720e+00
  NL step  1, |residual|_2 = 1.838827e-07
  NL step  2, |residual|_2 = 2.526653e-15
Writing state 92 at time 9.100000e-09 to file bih.1.e; output a total of 92 states so far
Solving for state 93, time 9.100000e-09
  NL step  0, |residual|_2 = 3.668992e+00
  NL step  1, |residual|_2 = 3.136350e-07
  NL step  2, |residual|_2 = 2.328935e-15
Writing state 93 at time 9.200000e-09 to file bih.1.e; output a total of 93 states so far
Solving for state 94, time 9.200000e-09
  NL step  0, |residual|_2 = 3.667269e+00
  NL step  1, |residual|_2 = 1.722002e-07
  NL step  2, |residual|_2 = 2.247037e-15
Writing state 94 at time 9.300000e-09 to file bih.1.e; output a total of 94 states so far
Solving for state 95, time 9.300000e-09
  NL step  0, |residual|_2 = 3.665550e+00
  NL step  1, |residual|_2 = 2.988151e-07
  NL step  2, |residual|_2 = 2.459460e-15
Writing state 95 at time 9.400000e-09 to file bih.1.e; output a total of 95 states so far
Solving for state 96, time 9.400000e-09
  NL step  0, |residual|_2 = 3.663837e+00
  NL step  1, |residual|_2 = 1.600913e-07
  NL step  2, |residual|_2 = 2.484313e-15
Writing state 96 at time 9.500000e-09 to file bih.1.e; output a total of 96 states so far
Solving for state 97, time 9.500000e-09
  NL step  0, |residual|_2 = 3.662129e+00
  NL step  1, |residual|_2 = 2.832423e-07
  NL step  2, |residual|_2 = 2.481569e-15
Writing state 97 at time 9.600000e-09 to file bih.1.e; output a total of 97 states so far
Solving for state 98, time 9.600000e-09
  NL step  0, |residual|_2 = 3.660426e+00
  NL step  1, |residual|_2 = 1.478468e-07
  NL step  2, |residual|_2 = 2.436396e-15
Writing state 98 at time 9.700000e-09 to file bih.1.e; output a total of 98 states so far
Solving for state 99, time 9.700000e-09
  NL step  0, |residual|_2 = 3.658728e+00
  NL step  1, |residual|_2 = 2.670648e-07
  NL step  2, |residual|_2 = 2.318336e-15
Writing state 99 at time 9.800000e-09 to file bih.1.e; output a total of 99 states so far
Solving for state 100, time 9.800000e-09
  NL step  0, |residual|_2 = 3.657035e+00
  NL step  1, |residual|_2 = 1.358089e-07
  NL step  2, |residual|_2 = 2.353539e-15
Writing state 100 at time 9.900000e-09 to file bih.1.e; output a total of 100 states so far
Solving for state 101, time 9.900000e-09
  NL step  0, |residual|_2 = 3.655347e+00
  NL step  1, |residual|_2 = 2.504283e-07
  NL step  2, |residual|_2 = 2.337847e-15
Writing state 101 at time 1.000000e-08 to file bih.1.e; output a total of 101 states so far
Solving for state 102, time 1.000000e-08
  NL step  0, |residual|_2 = 3.653664e+00
  NL step  1, |residual|_2 = 1.243929e-07
  NL step  2, |residual|_2 = 2.516733e-15
Writing state 102 at time 1.010000e-08 to file bih.1.e; output a total of 102 states so far
Writing state 102 at time 1.010000e-08 to file bih.1.e; output a total of 103 states so far

 -------------------------------------------------------------------------------------------------------------------
| Processor id:   0                                                                                                 |
| Num Processors: 4                                                                                                 |
| Time:           Fri Apr 19 11:51:51 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=12.2918, Active time=8.91372                                                   |
 ----------------------------------------------------------------------------------------------------------------
| 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.0010      0.001048    0.0011      0.001063    0.01     0.01     |
|   build_constraint_matrix()        120064    0.1420      0.000001    0.1420      0.000001    1.59     1.59     |
|   build_sparsity()                 1         0.0009      0.000905    0.0025      0.002509    0.01     0.03     |
|   constrain_elem_matrix()          47104     0.0514      0.000001    0.0514      0.000001    0.58     0.58     |
|   constrain_elem_vector()          72960     0.0755      0.000001    0.0755      0.000001    0.85     0.85     |
|   create_dof_constraints()         1         0.0121      0.012083    0.0191      0.019082    0.14     0.21     |
|   distribute_dofs()                1         0.0038      0.003789    0.0081      0.008060    0.04     0.09     |
|   dof_indices()                    146950    0.8103      0.000006    0.8103      0.000006    9.09     9.09     |
|   enforce_constraints_exactly()    470       0.1099      0.000234    0.1099      0.000234    1.23     1.23     |
|   prepare_send_list()              1         0.0000      0.000002    0.0000      0.000002    0.00     0.00     |
|   reinit()                         1         0.0039      0.003938    0.0039      0.003938    0.04     0.04     |
|                                                                                                                |
| EquationSystems                                                                                                |
|   build_solution_vector()          103       0.2317      0.002250    0.4340      0.004214    2.60     4.87     |
|                                                                                                                |
| ExodusII_IO                                                                                                    |
|   write_nodal_data()               103       5.3463      0.051906    5.3463      0.051906    59.98    59.98    |
|                                                                                                                |
| FE                                                                                                             |
|   compute_shape_functions()        473       0.0079      0.000017    0.0079      0.000017    0.09     0.09     |
|   init_shape_functions()           473       0.0165      0.000035    0.0165      0.000035    0.19     0.19     |
|   inverse_map()                    17        0.0001      0.000003    0.0001      0.000003    0.00     0.00     |
|                                                                                                                |
| FEMap                                                                                                          |
|   compute_affine_map()             120068    0.2075      0.000002    0.2075      0.000002    2.33     2.33     |
|   compute_face_map()               2         0.0000      0.000003    0.0000      0.000003    0.00     0.00     |
|   init_face_shape_functions()      2         0.0000      0.000003    0.0000      0.000003    0.00     0.00     |
|   init_reference_to_physical_map() 473       0.0023      0.000005    0.0023      0.000005    0.03     0.03     |
|                                                                                                                |
| Mesh                                                                                                           |
|   find_neighbors()                 1         0.0016      0.001646    0.0018      0.001840    0.02     0.02     |
|   renumber_nodes_and_elem()        2         0.0001      0.000062    0.0001      0.000062    0.00     0.00     |
|                                                                                                                |
| MeshCommunication                                                                                              |
|   compute_hilbert_indices()        2         0.0085      0.004269    0.0085      0.004269    0.10     0.10     |
|   find_global_indices()            2         0.0013      0.000633    0.0261      0.013068    0.01     0.29     |
|   parallel_sort()                  2         0.0045      0.002233    0.0144      0.007190    0.05     0.16     |
|                                                                                                                |
| MeshOutput                                                                                                     |
|   write_equation_systems()         103       0.0027      0.000026    5.7906      0.056219    0.03     64.96    |
|                                                                                                                |
| MeshTools::Generation                                                                                          |
|   build_cube()                     1         0.0006      0.000617    0.0006      0.000617    0.01     0.01     |
|                                                                                                                |
| MetisPartitioner                                                                                               |
|   partition()                      1         0.0035      0.003480    0.0176      0.017597    0.04     0.20     |
|                                                                                                                |
| Parallel                                                                                                       |
|   allgather()                      7         0.0010      0.000138    0.0013      0.000189    0.01     0.01     |
|   max(bool)                        1         0.0000      0.000004    0.0000      0.000004    0.00     0.00     |
|   max(scalar)                      5391      0.0257      0.000005    0.0257      0.000005    0.29     0.29     |
|   max(vector)                      1271      0.0158      0.000012    0.0336      0.000026    0.18     0.38     |
|   min(bool)                        6656      0.0297      0.000004    0.0297      0.000004    0.33     0.33     |
|   min(scalar)                      5389      0.0416      0.000008    0.0416      0.000008    0.47     0.47     |
|   min(vector)                      1271      0.0123      0.000010    0.0299      0.000023    0.14     0.34     |
|   probe()                          36        0.0014      0.000040    0.0014      0.000040    0.02     0.02     |
|   receive()                        36        0.0001      0.000003    0.0016      0.000044    0.00     0.02     |
|   send()                           36        0.0001      0.000002    0.0001      0.000002    0.00     0.00     |
|   send_receive()                   40        0.0002      0.000005    0.0019      0.000048    0.00     0.02     |
|   sum()                            785       0.0225      0.000029    0.0401      0.000051    0.25     0.45     |
|                                                                                                                |
| Parallel::Request                                                                                              |
|   wait()                           36        0.0001      0.000001    0.0001      0.000001    0.00     0.00     |
|                                                                                                                |
| Partitioner                                                                                                    |
|   set_node_processor_ids()         1         0.0007      0.000737    0.0063      0.006303    0.01     0.07     |
|   set_parent_processor_ids()       1         0.0003      0.000294    0.0003      0.000294    0.00     0.00     |
|                                                                                                                |
| PetscNonlinearSolver                                                                                           |
|   jacobian()                       184       0.6562      0.003566    1.2088      0.006569    7.36     13.56    |
|   residual()                       285       0.8731      0.003064    1.7258      0.006055    9.80     19.36    |
|   solve()                          101       0.1751      0.001734    3.1098      0.030790    1.96     34.89    |
|                                                                                                                |
| PointLocatorTree                                                                                               |
|   init(no master)                  1         0.0064      0.006400    0.0066      0.006645    0.07     0.07     |
|   operator()                       4         0.0001      0.000022    0.0001      0.000027    0.00     0.00     |
|                                                                                                                |
| System                                                                                                         |
|   project_vector()                 1         0.0033      0.003339    0.0053      0.005274    0.04     0.06     |
|   solve()                          101       0.0016      0.000015    3.1113      0.030805    0.02     34.91    |
 ----------------------------------------------------------------------------------------------------------------
| Totals:                            531016    8.9137                                          100.00            |
 ----------------------------------------------------------------------------------------------------------------

 
***************************************************************
* Done Running Example miscellaneous_ex7:
*  mpirun -np 4 example-devel --verbose dim=1 N=1024 initial_state=strip initial_center=0.5 initial_width=0.1 dt=1e-10 max_time=1e-8 -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_ex7'

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

Hosted By:
SourceForge.net Logo