def main(self, *args, **kwds):
    """
    Run the application.
    """
    if self.pdbOn:
          import pdb
          pdb.set_trace()
        
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()
    if 0 == comm.rank:
      self._info.log("Running on %d process(es)." % comm.size)

    from pylith.utils.profiling import resourceUsageString    
    self._debug.log(resourceUsageString())

    self._setupLogging()

    # Create mesh (adjust to account for interfaces (faults) if necessary)
    self._eventLogger.stagePush("Meshing")
    interfaces = None
    if "interfaces" in dir(self.problem):
      interfaces = self.problem.interfaces.components()
    mesh = self.mesher.create(self.problem.normalizer, interfaces)
    del interfaces
    del self.mesher
    self._debug.log(resourceUsageString())
    self._eventLogger.stagePop()

    # Setup problem, verify configuration, and then initialize
    self._eventLogger.stagePush("Setup")
    self.problem.preinitialize(mesh)
    self._debug.log(resourceUsageString())

    self.problem.verifyConfiguration()

    self.problem.initialize()
    self._debug.log(resourceUsageString())

    self._eventLogger.stagePop()

    # If initializing only, stop before running problem
    if self.initializeOnly:
      return

    # Run problem
    self.problem.run(self)
    self._debug.log(resourceUsageString())

    # Cleanup
    self._eventLogger.stagePush("Finalize")
    self.problem.finalize()
    self._eventLogger.stagePop()

    self.perfLogger.logMesh('Mesh', mesh)
    self.compilePerformanceLog()
    if self.perfLogger.verbose:
      self.perfLogger.show()

    return

  def _setupMaterials(self, materials):
    """
    Setup materials as integrators.
    """
    from pylith.feassemble.Integrator import implementsIntegrator
    
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    if 0 == comm.rank:
      self._info.log("Pre-initializing materials.")
    self._debug.log(resourceUsageString())
    for material in materials.components():
      integrator = self.elasticityIntegrator()
      if not implementsIntegrator(integrator):
        raise TypeError, \
              "Could not use '%s' as an integrator for material '%s'. " \
              "Functionality missing." % (integrator.name, material.label())
      integrator.preinitialize(self.mesh(), material)
      self.integrators.append(integrator)
      self._debug.log(resourceUsageString())

      if 0 == comm.rank:
        self._info.log("Added elasticity integrator for material '%s'." % material.label())
    return

  def _setupBC(self, boundaryConditions):
    """
    Setup boundary conditions as integrators or constraints.
    """
    from pylith.feassemble.Integrator import implementsIntegrator
    from pylith.feassemble.Constraint import implementsConstraint
    from pylith.bc.PointForce import PointForce

    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    if 0 == comm.rank:
      self._info.log("Pre-initializing boundary conditions.")
    self._debug.log(resourceUsageString())
    for bc in boundaryConditions.components():
      bc.preinitialize(self.mesh())
      foundType = False
      if implementsIntegrator(bc):
        foundType = True
        self.integrators.append(bc)
        if 0 == comm.rank:
          self._info.log("Added boundary condition '%s' as an integrator." % \
                           bc.label())
      if implementsConstraint(bc):
        foundType = True
        self.constraints.append(bc)
        if 0 == comm.rank:
          self._info.log("Added boundary condition '%s' as a constraint." % \
                           bc.label())
      if not foundType:
        raise TypeError, \
              "Could not determine whether boundary condition '%s' is an " \
              "integrator or a constraint." % bc.name
    self._debug.log(resourceUsageString())    
    return

  def finalize(self):
    """
    Cleanup after time stepping.
    """
    logEvent = "%sfinalize" % self._loggingPrefix
    self._eventLogger.eventBegin(logEvent)

    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    if 0 == comm.rank:
      self._info.log("Formulation finalize.")
    self._debug.log(resourceUsageString())
    for integrator in self.integrators:
      integrator.finalize()
    for constraint in self.constraints:
      constraint.finalize()
    for output in self.output.components():
      output.close()
      output.finalize()
    self._debug.log(resourceUsageString())
    
    self._modelMemoryUse()

    self._eventLogger.eventEnd(logEvent)
    return

  def create(self, normalizer, faults=None):
    """
    Hook for creating mesh.
    """
    from pylith.utils.profiling import resourceUsageString
    from pylith.mpi.Communicator import petsc_comm_world
    comm = petsc_comm_world()

    self._setupLogging()
    logEvent = "%screate" % self._loggingPrefix
    self._eventLogger.eventBegin(logEvent)    

    # Read mesh
    mesh = self.reader.read(self.debug, self.interpolate)
    if self.debug:
      mesh.view()

    # Reorder mesh
    if self.reorderMesh:
      logEvent2 = "%sreorder" % self._loggingPrefix
      self._eventLogger.eventBegin(logEvent2)
      self._debug.log(resourceUsageString())
      if 0 == comm.rank:
        self._info.log("Reordering cells and vertices.")
      from pylith.topology.ReverseCuthillMcKee import ReverseCuthillMcKee
      ordering = ReverseCuthillMcKee()
      ordering.reorder(mesh)
      self._eventLogger.eventEnd(logEvent2)

    # Adjust topology
    self._debug.log(resourceUsageString())
    if 0 == comm.rank:
      self._info.log("Adjusting topology.")
    self._adjustTopology(mesh, faults)

    # Distribute mesh
    if comm.size > 1:
      if 0 == comm.rank:
        self._info.log("Distributing mesh.")
      mesh = self.distributor.distribute(mesh, normalizer)
      if self.debug:
        mesh.view()
      mesh.memLoggingStage = "DistributedMesh"

    # Refine mesh (if necessary)
    newMesh = self.refiner.refine(mesh)
    if not newMesh == mesh:
      mesh.cleanup()
      newMesh.memLoggingStage = "RefinedMesh"

    # Can't reorder mesh again, because we do not have routine to
    # unmix normal and hybrid cells.

    # Nondimensionalize mesh (coordinates of vertices).
    from pylith.topology.topology import MeshOps_nondimensionalize
    MeshOps_nondimensionalize(newMesh, normalizer)

    self._eventLogger.eventEnd(logEvent)    
    return newMesh

  def initialize(self, dimension, normalizer):
    """
    Initialize problem for implicit time integration.
    """
    logEvent = "%sinit" % self._loggingPrefix
    self._eventLogger.eventBegin(logEvent)

    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    self._initialize(dimension, normalizer)

    #from pylith.utils.petsc import MemoryLogger
    #memoryLogger = MemoryLogger.singleton()
    #memoryLogger.setDebug(0)
    #memoryLogger.stagePush("Problem")

    # Allocate other fields, reusing layout from dispIncr
    if 0 == comm.rank:
      self._info.log("Creating other fields.")
    self.fields.add("velocity(t)", "velocity")
    self.fields.copyLayout("dispIncr(t->t+dt)")

    # Setup fields and set to zero
    dispT = self.fields.get("disp(t)")
    dispT.zeroAll()
    residual = self.fields.get("residual")
    residual.zeroAll()
    residual.createScatter(residual.mesh())

    lengthScale = normalizer.lengthScale()
    timeScale = normalizer.timeScale()
    velocityScale = lengthScale / timeScale
    velocityT = self.fields.get("velocity(t)")
    velocityT.scale(velocityScale.value)
    velocityT.zeroAll()

    self._debug.log(resourceUsageString())
    #memoryLogger.stagePop()

    # Allocates memory for nonzero pattern and Jacobian
    if 0 == comm.rank:
      self._info.log("Creating Jacobian matrix.")
    self._setJacobianMatrixType()
    from pylith.topology.Jacobian import Jacobian
    self.jacobian = Jacobian(self.fields.solution(),
                             self.matrixType, self.blockMatrixOkay)
    self.jacobian.zero() # TEMPORARY, to get correct memory usage
    self._debug.log(resourceUsageString())

    #memoryLogger.stagePush("Problem")
    if 0 == comm.rank:
      self._info.log("Initializing solver.")
    self.solver.initialize(self.fields, self.jacobian, self)
    self._debug.log(resourceUsageString())

    #memoryLogger.stagePop()
    #memoryLogger.setDebug(0)
    return

  def _setupInterfaces(self, interfaceConditions):
    """
    Setup interfaces as integrators or constraints.
    """
    from pylith.feassemble.Integrator import implementsIntegrator
    from pylith.feassemble.Constraint import implementsConstraint

    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    if 0 == comm.rank:
      self._info.log("Pre-initializing interior interfaces.")
    for ic in interfaceConditions.components():
      ic.preinitialize(self.mesh())
      foundType = False
      if implementsIntegrator(ic):
        foundType = True
        self.integrators.append(ic)
        if 0 == comm.rank:
          self._info.log("Added interface condition '%s' as an integrator." % \
                           ic.label())
      if implementsConstraint(ic):
        foundType = True
        self.constraints.append(ic)
        if 0 == comm.rank:
          self._info.log("Added interface condition '%s' as a constraint." % \
                           ic.label())
      if not foundType:
        raise TypeError, \
              "Could not determine whether interface condition '%s' is an " \
              "integrator or a constraint." % ic.name
    self._debug.log(resourceUsageString())    
    return

  def main(self, *args, **kwds):
    """
    Run the application.
    """
    from pylith.utils.profiling import resourceUsageString
    
    self.petsc.initialize()
    self._debug.log(resourceUsageString())

    # Create mesh (adjust to account for interfaces (faults) if necessary)
    interfaces = None
    mesh = self.mesher.create(interfaces)
    self._debug.log(resourceUsageString())

    self.petsc.finalize()
    return

  def _reformJacobian(self, t, dt):
    """
    Reform Jacobian matrix for operator.
    """
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    self._debug.log(resourceUsageString())
    if 0 == comm.rank:
      self._info.log("Integrating Jacobian operator.")
    self._eventLogger.stagePush("Reform Jacobian")

    self.updateSettings(self.jacobian, self.fields, t, dt)
    ModuleExplicit.reformJacobianLumped(self)

    self._eventLogger.stagePop()

    if self.viewJacobian:
      self.jacobian.view("Jacobian")

    self._debug.log(resourceUsageString())
    return

  def _reformJacobian(self, t, dt):
    """
    Reform Jacobian matrix for operator.
    """
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    self._debug.log(resourceUsageString())
    if 0 == comm.rank:
      self._info.log("Integrating Jacobian operator.")
    self._eventLogger.stagePush("Reform Jacobian")

    self.updateSettings(self.jacobian, self.fields, t, dt)
    self.reformJacobian()

    self._eventLogger.stagePop()

    if self.viewJacobian:
      from pylith.mpi.Communicator import Communicator
      comm = Communicator(self.mesh().comm())
      self.jacobianViewer.view(self.jacobian, t, comm)

    self._debug.log(resourceUsageString())
    return

  def _reformResidual(self, t, dt):
    """
    Reform residual vector for operator.
    """
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    if 0 == comm.rank:
      self._info.log("Integrating residual term in operator.")
    self._eventLogger.stagePush("Reform Residual")

    self.updateSettings(self.jacobian, self.fields, t, dt)
    self.reformResidual()

    self._eventLogger.stagePop()
    self._debug.log(resourceUsageString())
    return

  def _showStatus(self, stage):
    import sys
    from pylith.utils.profiling import resourceUsageString
    from pylith.mpi.Communicator import petsc_comm_world
    comm = petsc_comm_world()

    if comm.rank == 0:
      print "\n----------------------------------------------------------------------"
      print "STAGE: %s" % stage
      print "----------------------------------------------------------------------"
    for irank in xrange(comm.size):
      comm.barrier()
      if comm.rank == irank:
        print "\nPROCESSOR %d" % comm.rank
        print "\nStatus from ps: %s\n" % resourceUsageString()
        self.logger.show()
        sys.stdout.flush()

      comm.barrier()
    return

  def create(self, normalizer, faults=None):
    """
    Hook for creating mesh.
    """
    from pylith.utils.profiling import resourceUsageString

    self._setupLogging()
    logEvent = "%screate" % self._loggingPrefix
    self._eventLogger.eventBegin(logEvent)    

    mesh = self.reader.read(self.debug, self.interpolate)
    if self.debug:
      mesh.view("Finite-element mesh.")
    self._debug.log(resourceUsageString())

    # refine mesh (if necessary)
    mesh = self.refiner.refine(mesh)

    # Nondimensionalize mesh (coordinates of vertices).
    from pylith.topology.topology import MeshOps_nondimensionalize
    MeshOps_nondimensionalize(mesh, normalizer)

    self._eventLogger.eventEnd(logEvent)    
    return mesh

  def initialize(self, dimension, normalizer):
    """
    Initialize problem for explicit time integration.
    """
    logEvent = "%sinit" % self._loggingPrefix
    self._eventLogger.eventBegin(logEvent)
    
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    self._initialize(dimension, normalizer)

    #from pylith.utils.petsc import MemoryLogger
    #memoryLogger = MemoryLogger.singleton()
    #memoryLogger.setDebug(0)
    #memoryLogger.stagePush("Problem")

    # Allocate other fields, reusing layout from dispIncr
    if 0 == comm.rank:
      self._info.log("Creating other fields.")
    self.fields.add("disp(t-dt)", "displacement")
    self.fields.add("velocity(t)", "velocity")
    self.fields.add("acceleration(t)", "acceleration")
    self.fields.copyLayout("dispIncr(t->t+dt)")
    self._debug.log(resourceUsageString())

    # Setup fields and set to zero
    dispTmdt = self.fields.get("disp(t-dt)")
    dispTmdt.zeroAll()
    dispT = self.fields.get("disp(t)")
    dispT.zeroAll()
    residual = self.fields.get("residual")
    residual.zeroAll()
    residual.createScatter(residual.mesh())

    lengthScale = normalizer.lengthScale()
    timeScale = normalizer.timeScale()
    velocityScale = lengthScale / timeScale
    velocityT = self.fields.get("velocity(t)")
    velocityT.scale(velocityScale.value)
    velocityT.zeroAll()

    accelerationScale = velocityScale / timeScale
    accelerationT = self.fields.get("acceleration(t)")
    accelerationT.scale(accelerationScale.value)
    accelerationT.zeroAll()

    self._debug.log(resourceUsageString())
    #memoryLogger.stagePop()

    if 0 == comm.rank:
      self._info.log("Creating lumped Jacobian matrix.")
    from pylith.topology.Field import Field
    jacobian = Field(self.mesh())
    jacobian.label("jacobian")

    # Setup section manually. Cloning the solution field includes
    # constraints which messes up the solve for constrained DOF.
    pressureScale = normalizer.pressureScale()
    jacobian.subfieldAdd("displacement", dimension, jacobian.VECTOR, lengthScale.value)
    jacobian.subfieldAdd("lagrange_multiplier", dimension, jacobian.VECTOR, pressureScale.value)
    jacobian.subfieldsSetup()
    jacobian.setupSolnChart()
    jacobian.setupSolnDof(dimension)
    # Loop over integrators to adjust DOF layout
    for integrator in self.integrators:
      integrator.setupSolnDof(jacobian)
    jacobian.vectorFieldType(jacobian.VECTOR)
    jacobian.allocate()
    jacobian.zeroAll()
    self.jacobian = jacobian
    self._debug.log(resourceUsageString())

    #memoryLogger.stagePush("Problem")
    if 0 == comm.rank:
      self._info.log("Initializing solver.")
    self.solver.initialize(self.fields, self.jacobian, self)
    self._debug.log(resourceUsageString())

    #memoryLogger.stagePop()
    #memoryLogger.setDebug(0)
    self._eventLogger.eventEnd(logEvent)
    return

  def _initialize(self, dimension, normalizer):
    """
    Create integrators for each element family.
    """
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    self.timeStep.initialize(normalizer)

    numTimeSteps = self.timeStep.numTimeSteps()
    totalTime = self.timeStep.totalTime

    from pylith.topology.SolutionFields import SolutionFields
    self.fields = SolutionFields(self.mesh())
    self._debug.log(resourceUsageString())

    if 0 == comm.rank:
      self._info.log("Initializing integrators.")
    for integrator in self.integrators:
      if not self.gravityField is None:
        integrator.gravityField(self.gravityField)
      integrator.initialize(totalTime, numTimeSteps, normalizer)
    ModuleFormulation.integrators(self, self.integrators)
    self._debug.log(resourceUsageString())

    if 0 == comm.rank:
      self._info.log("Initializing constraints.")
    for constraint in self.constraints:
      constraint.initialize(totalTime, numTimeSteps, normalizer)
    self._debug.log(resourceUsageString())

    if 0 == comm.rank:
      self._info.log("Setting up solution output.")
    for output in self.output.components():
      output.initialize(self.mesh(), normalizer)
      output.writeInfo()
      output.open(totalTime, numTimeSteps)
    self._debug.log(resourceUsageString())

    # Setup fields
    if 0 == comm.rank:
      self._info.log("Creating solution field.")
    #from pylith.utils.petsc import MemoryLogger
    #memoryLogger = MemoryLogger.singleton()
    #memoryLogger.setDebug(0)
    #memoryLogger.stagePush("Problem")
    self.fields.add("dispIncr(t->t+dt)", "displacement_increment")
    self.fields.add("disp(t)", "displacement")
    self.fields.add("residual", "residual")
    self.fields.solutionName("dispIncr(t->t+dt)")

    lengthScale = normalizer.lengthScale()
    pressureScale = normalizer.pressureScale()

    solution = self.fields.get("dispIncr(t->t+dt)")
    solution.subfieldAdd("displacement", dimension, solution.VECTOR, lengthScale.value)
    solution.subfieldAdd("lagrange_multiplier", dimension, solution.VECTOR, pressureScale.value)
    solution.subfieldsSetup()
    solution.setupSolnChart()
    solution.setupSolnDof(dimension)
    # Loop over integrators to adjust DOF layout
    for integrator in self.integrators:
      integrator.setupSolnDof(solution)
    solution.vectorFieldType(solution.VECTOR)
    solution.scale(lengthScale.value)

    for constraint in self.constraints:
      constraint.setConstraintSizes(solution)
    solution.allocate()
    solution.zeroAll()
    for constraint in self.constraints:
      constraint.setConstraints(solution)
    for integrator in self.integrators:
      integrator.checkConstraints(solution)

    #memoryLogger.stagePop()

    # This also creates a global order.
    solution.createScatter(solution.mesh())

    #memoryLogger.stagePush("Problem")
    dispT = self.fields.get("disp(t)")
    dispT.vectorFieldType(dispT.VECTOR)
    dispT.scale(lengthScale.value)

    residual = self.fields.get("residual")
    residual.vectorFieldType(residual.VECTOR)
    residual.scale(lengthScale.value)

    #memoryLogger.stagePop()
    #memoryLogger.setDebug(0)
    self._debug.log(resourceUsageString())

    return

  def initialize(self, dimension, normalizer):
    """
    Initialize problem for explicit time integration.
    """
    logEvent = "%sinit" % self._loggingPrefix
    self._eventLogger.eventBegin(logEvent)
    
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    self._initialize(dimension, normalizer)

    #from pylith.utils.petsc import MemoryLogger
    #memoryLogger = MemoryLogger.singleton()
    #memoryLogger.setDebug(0)
    #memoryLogger.stagePush("Problem")

    # Allocate other fields, reusing layout from dispIncr
    if 0 == comm.rank:
      self._info.log("Creating other fields.")
    self.fields.add("disp(t-dt)", "displacement")
    self.fields.add("velocity(t)", "velocity")
    self.fields.add("acceleration(t)", "acceleration")
    self.fields.copyLayout("dispIncr(t->t+dt)")
    self._debug.log(resourceUsageString())

    # Setup fields and set to zero
    dispTmdt = self.fields.get("disp(t-dt)")
    dispTmdt.zeroAll()
    dispT = self.fields.get("disp(t)")
    dispT.zeroAll()
    residual = self.fields.get("residual")
    residual.zeroAll()
    residual.createScatter(residual.mesh())

    lengthScale = normalizer.lengthScale()
    timeScale = normalizer.timeScale()
    velocityScale = lengthScale / timeScale
    velocityT = self.fields.get("velocity(t)")
    velocityT.scale(velocityScale.value)
    velocityT.zeroAll()

    accelerationScale = velocityScale / timeScale
    accelerationT = self.fields.get("acceleration(t)")
    accelerationT.scale(accelerationScale.value)
    accelerationT.zeroAll()

    self._debug.log(resourceUsageString())
    #memoryLogger.stagePop()

    if 0 == comm.rank:
      self._info.log("Creating lumped Jacobian matrix.")
    from pylith.topology.Field import Field
    jacobian = Field(self.mesh())
    jacobian.newSection(jacobian.VERTICES_FIELD, dimension)
    jacobian.allocate()
    jacobian.label("jacobian")
    jacobian.vectorFieldType(jacobian.VECTOR)
    self.jacobian = jacobian
    self._debug.log(resourceUsageString())

    #memoryLogger.stagePush("Problem")
    if 0 == comm.rank:
      self._info.log("Initializing solver.")
    self.solver.initialize(self.fields, self.jacobian, self)
    self._debug.log(resourceUsageString())

    #memoryLogger.stagePop()
    #memoryLogger.setDebug(0)
    self._eventLogger.eventEnd(logEvent)
    return