def preinitialize(self, mesh):
    """
    Do pre-initialization setup.
    """
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()    

    if 0 == comm.rank:
      self._info.log("Pre-initializing fault '%s'." % self.label())
    FaultCohesive.preinitialize(self, mesh)
    Integrator.preinitialize(self, mesh)

    ModuleFaultCohesiveDyn.quadrature(self, self.faultQuadrature)

    if mesh.dimension() == 2:
      self.availableFields['vertex']['info'] += ["strike_dir"]
    elif mesh.dimension() == 3:
      self.availableFields['vertex']['info'] += ["strike_dir",
                                                 "dip_dir"]

    if not isinstance(self.tract, NullComponent):
      self.tract.preinitialize(mesh)
      self.availableFields['vertex']['info'] += self.tract.availableFields['vertex']['info']

    self.availableFields['vertex']['info'] += \
        self.friction.availableFields['vertex']['info']
    self.availableFields['vertex']['data'] += \
        self.friction.availableFields['vertex']['data']
    return

  def preinitialize(self, mesh):
    """
    Setup integrators for each element family (material/quadrature,
    bc/quadrature, etc.).
    """
    self._setupLogging()
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()
    
    if 0 == comm.rank:
      self._info.log("Pre-initializing problem.")
    import weakref
    self.mesh = weakref.ref(mesh)
    self.formulation.preinitialize(mesh, self.materials, self.bc, self.interfaces, self.gravityField)

    # Find fault for impulses
    found = False
    for fault in self.interfaces.components():
      if self.faultId == fault.id():
        self.source = fault
        found = True
        break
    if not found:
      raise ValueError("Could not find fault interface with id '%d' for "
                       "Green's function impulses." % self.faultId)
    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 read(self, debug, interpolate):
    """
    Read finite-element mesh and store in Sieve mesh object.

    @returns PETSc mesh object containing finite-element mesh
    """
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()
    if 0 == comm.rank:
      self._info.log("Reading finite-element mesh")

    # Set flags
    self.debug(debug)
    self.interpolate(interpolate)

    # Initialize coordinate system
    if self.coordsys is None:
      raise ValueError, "Coordinate system for mesh is unknown."

    from pylith.mpi.Communicator import petsc_comm_world
    from pylith.topology.Mesh import Mesh    
    mesh = Mesh(dim=self.coordsys.spaceDim(), comm=petsc_comm_world())
    mesh.coordsys(self.coordsys)

    # Read mesh
    ModuleMeshIO.read(self, mesh)
    return mesh

  def preinitialize(self, mesh):
    """
    Do pre-initialization setup.
    """
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    if 0 == comm.rank:
      self._info.log("Pre-initializing fault '%s'." % self.label())
    FaultCohesive.preinitialize(self, mesh)
    Integrator.preinitialize(self, mesh)

    ModuleFaultCohesiveKin.quadrature(self, self.faultQuadrature)

    for eqsrc in self.eqsrcs.components():
      eqsrc.preinitialize()
    ModuleFaultCohesiveKin.eqsrcs(self, self.eqsrcs.inventory.facilityNames(),
                                  self.eqsrcs.components())

    for name in self.eqsrcs.inventory.facilityNames():
      self.availableFields['vertex']['info'] += ["final_slip_%s" % name]
      self.availableFields['vertex']['info'] += ["slip_time_%s" % name]

    if mesh.dimension() == 2:
      self.availableFields['vertex']['info'] += ["strike_dir"]
    elif mesh.dimension() == 3:
      self.availableFields['vertex']['info'] += ["strike_dir",
                                                 "dip_dir"]
    return

  def run(self, app):
    """
    Compute Green's functions associated with fault slip.
    """
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    if 0 == comm.rank:
      self._info.log("Computing Green's functions.")
    self.checkpointTimer.toplevel = app # Set handle for saving state

    # Limit material behavior to linear regime
    for material in self.materials.components():
      material.useElasticBehavior(True)

    nimpulses = self.source.numImpulses()
    if nimpulses > 0:
      self.progressMonitor.open()
    
    ipulse = 0;
    dt = 1.0
    while ipulse < nimpulses:
      self.progressMonitor.update(ipulse, 0, nimpulses)

      self._eventLogger.stagePush("Prestep")
      if 0 == comm.rank:
        self._info.log("Main loop, impulse %d of %d." % (ipulse+1, nimpulses))
      
      # Implicit time stepping computes solution at t+dt, so set
      # t=ipulse-dt, so that t+dt corresponds to the impulse
      t = float(ipulse)-dt

      # Checkpoint if necessary
      self.checkpointTimer.update(t)

      if 0 == comm.rank:
        self._info.log("Preparing impulse %d of %d." % \
                         (ipulse+1, nimpulses))
      self.formulation.prestep(t, dt)
      self._eventLogger.stagePop()

      if 0 == comm.rank:
        self._info.log("Computing response to impulse %d of %d." %
                         (ipulse+1, nimpulses))
      self._eventLogger.stagePush("Step")
      self.formulation.step(t, dt)
      self._eventLogger.stagePop()

      if 0 == comm.rank:
        self._info.log("Finishing impulse %d of %d." % \
                         (ipulse+1, nimpulses))
      self._eventLogger.stagePush("Poststep")
      self.formulation.poststep(t, dt)
      self._eventLogger.stagePop()

      # Update time/impulse
      ipulse += 1

    self.progressMonitor.close()      
    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.
    """
    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 poststep(self, t, dt):
    """
    Hook for doing stuff after advancing time step.
    """
    logEvent = "%spoststep" % self._loggingPrefix
    self._eventLogger.eventBegin(logEvent)
    
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    # The velocity and acceleration at time t depends on the
    # displacement at time t+dt, we want to output BEFORE updating the
    # displacement fields so that the displacement, velocity, and
    # acceleration files are all at time t.
    if 0 == comm.rank:
      self._info.log("Writing solution fields.")
    for output in self.output.components():
      output.writeData(t, self.fields)
    self._writeData(t)

    # Update displacement field from time t to time t+dt.
    dispIncr = self.fields.get("dispIncr(t->t+dt)")
    dispT = self.fields.get("disp(t)")
    dispTmdt = self.fields.get("disp(t-dt)")

    dispTmdt.copy(dispT)
    dispT.add(dispIncr)
    dispIncr.zeroAll()

    # Complete post-step processing.
    Formulation.poststep(self, t, dt)

    self._eventLogger.eventEnd(logEvent)    
    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 preinitialize(self, mesh, materials, boundaryConditions,
                    interfaceConditions, gravityField):
    """
    Create integrator for each element family.
    """
    self._setupLogging()
    logEvent = "%spreinit" % self._loggingPrefix
    self._eventLogger.eventBegin(logEvent)

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

    self.timeStep.preinitialize()

    import weakref
    self.mesh = weakref.ref(mesh)
    self.integrators = []
    self.constraints = []
    self.gravityField = gravityField

    self.solver.preinitialize()
    self._setupMaterials(materials)
    self._setupBC(boundaryConditions)
    self._setupInterfaces(interfaceConditions)

    if 0 == comm.rank:
      self._info.log("Pre-initializing output.")
    for output in self.output.components():
      output.preinitialize()

    self._eventLogger.eventEnd(logEvent)
    return

  def poststep(self, t, dt):
    """
    Hook for doing stuff after advancing time step.
    """
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    # Update displacement field from time t to time t+dt.
    dispIncr = self.fields.get("dispIncr(t->t+dt)")
    disp = self.fields.get("disp(t)")
    disp.add(dispIncr)
    dispIncr.zeroAll()

    # Complete post-step processing, then write data.
    Formulation.poststep(self, t, dt)

    # Write data. Velocity at time t will be based upon displacement
    # at time t-dt and t.
    if 0 == comm.rank:
      self._info.log("Writing solution fields.")
    for output in self.output.components():
      output.writeData(t+dt, self.fields)
    self._writeData(t+dt)

    return

  def step(self, t, dt):
    """
    Advance to next time step.
    """
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    dispIncr = self.fields.get("dispIncr(t->t+dt)")

    self._reformResidual(t+dt, dt)

    if 0 == comm.rank:
      self._info.log("Solving equations.")
    self._eventLogger.stagePush("Solve")

    residual = self.fields.get("residual")
    #self.jacobian.view() # TEMPORARY
    self.solver.solve(dispIncr, self.jacobian, residual)
    #dispIncr.view("DISP INCR") # TEMPORARY

    # DEBUGGING Verify solution makes residual 0
    #self._reformResidual(t+dt, dt)
    #residual.view("RESIDUAL")
    
    self._eventLogger.stagePop()

    return

  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 resourceUsageString():
  """
  Get CPU time and memory usage as a string.
  """
  from pylith.mpi.Communicator import mpi_comm_world
  comm = mpi_comm_world()
  (cputime, memory) = resourceUsage()
  return "[%d] CPU time: %s, Memory usage: %.2f MB" % \
      (comm.rank, cputime, memory)

  def write(self, mesh):
    """
    Write finite-element mesh.stored in Sieve mesh object.

    @param mesh PETSc mesh object containing finite-element mesh
    """
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()
    if 0 == comm.rank:
      self._info.log("Writing finite-element mesh")
    ModuleMeshIO.write(self, mesh)
    return

 def _createPath(self, filename):
     """Create path for filename if it doesn't exist.
     """
     import os
     relpath = os.path.dirname(filename)
 
     if len(relpath) > 0 and not os.path.exists(relpath):
         # Only create directory on proc 0
         from pylith.mpi.Communicator import mpi_comm_world
         comm = mpi_comm_world()
         if 0 == comm.rank:
             os.makedirs(relpath)
     return

  def initialize(self):
    """
    Setup integrators for each element family (material/quadrature,
    bc/quadrature, etc.).
    """
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()

    if 0 == comm.rank:
      self._info.log("Initializing problem.")
    self.checkpointTimer.initialize(self.normalizer)
    self.formulation.initialize(self.dimension, self.normalizer)
    return

  def close(self):
    """
    Close writer.
    """
    ModuleDataWriterHDF5Ext.close(self)

    # Only write Xdmf file on proc 0
    from pylith.mpi.Communicator import mpi_comm_world
    comm = mpi_comm_world()
    if not comm.rank:
      from Xdmf import Xdmf
      xdmf = Xdmf()
      xdmf.write(ModuleDataWriterHDF5Ext.hdf5Filename(self), verbose=False)
    return