def setUp(self):
    """
    Setup mesh and associated field.
    """
    from spatialdata.geocoords.CSCart import CSCart
    cs = CSCart()
    cs.inventory.spaceDim = 2
    cs._configure()

    from spatialdata.units.Nondimensional import Nondimensional
    normalizer = Nondimensional()
    normalizer._configure()    

    from pylith.meshio.MeshIOAscii import MeshIOAscii
    importer = MeshIOAscii()
    importer.inventory.filename = "data/tri3.mesh"
    importer.inventory.coordsys = cs
    importer._configure()
    self.mesh = importer.read(debug=False, interpolate=False)

    from pylith.topology.SolutionFields import SolutionFields
    fields = SolutionFields(self.mesh)
    fields.add("disp t+dt", "displacement")
    fields.solutionName("disp t+dt")
    solution = fields.solution()
    solution.newSection(solution.VERTICES_FIELD, self.mesh.dimension())
    solution.allocate()
    solution.zero()

    self.fields = fields
    self.jacobian = Jacobian(solution)
    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 test_integrateJacobian(self):
    """
    Test integrateJacobian().

    This does nothing for Neumann BC.
    """

    (mesh, bc, fields) = self._initialize()

    from pylith.topology.Jacobian import Jacobian
    jacobian = Jacobian(fields.solution())
    jacobian.zero()
    t = 0.24
    bc.integrateJacobian(jacobian, t, fields)
    self.assertEqual(False, bc.needNewJacobian())

    # No testing of result.
    return

  def test_integrateJacobian(self):
    """
    Test integrateJacobian().

    WARNING: This is not a rigorous test of integrateJacobian() because we
    don't verify the results.
    """

    (mesh, bc, fields) = self._initialize()

    from pylith.topology.Jacobian import Jacobian
    jacobian = Jacobian(fields.solution())
    jacobian.zero()
    t = 0.24
    bc.integrateJacobian(jacobian, t, fields)
    self.assertEqual(False, bc.needNewJacobian())

    # No testing of result.
    return

  def test_integrateJacobian(self):
    """
    Test integrateJacobian().

    WARNING: This is not a rigorous test of integrateJacobian() because we
    neither set the input fields or verify the results.
    """
    (mesh, fault, fields) = self._initialize()

    from pylith.topology.Jacobian import Jacobian
    jacobian = Jacobian(fields.solution())
    jacobian.zero()
    t = 1.0
    fault.integrateJacobian(jacobian, t, fields)
    self.assertEqual(False, fault.needNewJacobian())

    # We should really add something here to check to make sure things
    # actually initialized correctly    
    return

  def test_integrateJacobian(self):
    """
    Test integrateJacobian().

    WARNING: This is not a rigorous test of integrateJacobian()
    because we neither set the input fields or verify the results.
    """
    (mesh, integrator) = self._preinitialize()
    fields = self._initialize(mesh, integrator)

    from pylith.topology.Jacobian import Jacobian
    jacobian = Jacobian(fields.solution())
    jacobian.zero()
    t = 7.3
    self.assertEqual(True, integrator.needNewJacobian())
    integrator.integrateJacobian(jacobian, t, fields)
    self.assertEqual(False, integrator.needNewJacobian())
    
    # No test of result.
    return

  def test_write(self):
    """
    Test write().

    :WARNING: This is not a complete test of write(). We do not
    verify the results.
    """
    self.jacobian = Jacobian(self.fields.solution(), "aij")
    self.jacobian.assemble("final_assembly")

    self.jacobian.write("jacobian.mat", self.mesh.comm())

    # No testing of result.
    return

class TestJacobian(unittest.TestCase):
  """
  Unit testing of Jacobian object.
  """

  def setUp(self):
    """
    Setup mesh and associated field.
    """
    from spatialdata.geocoords.CSCart import CSCart
    cs = CSCart()
    cs.inventory.spaceDim = 2
    cs._configure()

    from spatialdata.units.Nondimensional import Nondimensional
    normalizer = Nondimensional()
    normalizer._configure()    

    from pylith.meshio.MeshIOAscii import MeshIOAscii
    importer = MeshIOAscii()
    importer.inventory.filename = "data/tri3.mesh"
    importer.inventory.coordsys = cs
    importer._configure()
    self.mesh = importer.read(debug=False, interpolate=False)

    from pylith.topology.SolutionFields import SolutionFields
    fields = SolutionFields(self.mesh)
    fields.add("disp t+dt", "displacement")
    fields.solutionName("disp t+dt")
    solution = fields.solution()
    solution.newSection(solution.VERTICES_FIELD, self.mesh.dimension())
    solution.allocate()
    solution.zero()

    self.fields = fields
    self.jacobian = Jacobian(solution)
    return


  def tearDown(self):
    self.jacobian.cleanup()
    return


  def test_constructor(self):
    """
    Test constructor.
    """
    # setUp() tests constructor with default type
    jacobianA = Jacobian(self.fields.solution(), "aij")
    jacobianB = Jacobian(self.fields.solution(), "baij")
    return


  def test_matrix(self):
    """
    Test matrix().

    :WARNING: This is not a complete test of matrix(). We do not
    verify the results.
    """
    matrix = self.jacobian.matrix()

    # No testing of result.
    return


  def test_assemble(self):
    """
    Test assemble().

    :WARNING: This is not a complete test of assemble(). We do not
    verify the results.
    """
    self.jacobian.assemble("flush_assembly")
    self.jacobian.assemble("final_assembly")

    # No testing of result.
    return


  def test_zero(self):
    """
    Test zero().

    :WARNING: This is not a complete test of zero(). We do not
    verify the results.
    """
    self.jacobian.zero()

    # No testing of result.
    return


  def test_view(self):
    """
    Test view().

    :WARNING: This is not a complete test of view(). We do not
    verify the results.
#.........这里部分代码省略.........

class Implicit(Formulation, ModuleImplicit):
  """
  Python Implicit object for solving equations using an implicit
  formulation.

  The formulation has the general form,

  [A(t+dt)] {u(t+dt)} = {b(t+dt)}. 

  We know the solution at time t, so we write {u(t+dt)} as {u(t)} +
  {du(t)}, where {du(t)} is the increment in the solution from time t
  to time t+dt. Thus, we solve

  [A(t+dt)] {du(t)} = {b(t+dt)} - [A(t+dt)]{u(t)}.

  We solve this system by forming the Jacobian, A, and the residual

  {r(t+dt)} = {b(t+dt)} - [A(t+dt)]{u(t)} - [A(t+dt)]{du(t)}

  which we combine into

  {r(t+dt)} = {b(t+dt)} - [A(t+dt)]{u(t)+du(t)}.

  The method reformJacobian() computes [A(t+dt)] and the method
  reformResidual computes {r(t+dt)}. Note that in forming the residual
  we compute the action [A(t+dt)]{u(t)+du(t)} and do not perform a
  matrix-vector multiplication.

  [A(t+dt)] generally depends on {u(t+dt)} as well as the current
  stresses and additional state variables.  

  For linear elastic or viscoelastic problems with constant time step
  size, A is a constant (after the elastic solution).  {b(t+dt)}
  generally depends on the loads applied for time step t+dt (including
  the contributions to the internal force vector from
  displacement/velocity BC) as well as the internal force vector
  computed from the current stresses.

  Factory: pde_formulation.
  """

  # INVENTORY //////////////////////////////////////////////////////////

  class Inventory(Formulation.Inventory):
    """
    Python object for managing Implicit facilities and properties.
    """

    ## @class Inventory
    ## Python object for managing Implicit facilities and properties.
    ##
    ## \b Properties
    ## @li None
    ##
    ## \b Facilities
    ## @li None

    import pyre.inventory


  # PUBLIC METHODS /////////////////////////////////////////////////////

  def __init__(self, name="implicit"):
    """
    Constructor.
    """
    Formulation.__init__(self, name)
    ModuleImplicit.__init__(self)
    self._loggingPrefix = "TSIm "
    return


  def elasticityIntegrator(self):
    """
    Get integrator for elastic material.
    """
    from pylith.feassemble.ElasticityImplicit import ElasticityImplicit
    integrator = ElasticityImplicit()
    return integrator


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

    comm = self.mesh().comm()

    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.")
#.........这里部分代码省略.........