def test_preinitialize(self):
    """
    Test preinitialize().

    WARNING: This is not a rigorous test of initialize() because we
    don't verify the results.
    """
    from pylith.feassemble.FIATSimplex import FIATSimplex
    cell = FIATSimplex()
    cell.inventory.dimension = 2
    cell.inventory.order = 1
    cell.inventory.degree = 1
    cell._configure()

    from pylith.feassemble.Quadrature import Quadrature
    quadrature = Quadrature()
    quadrature.inventory.cell = cell
    quadrature.inventory.minJacobian = 1.0e-4
    quadrature._configure()

    from spatialdata.spatialdb.SimpleDB import SimpleDB
    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
    iohandler = SimpleIOAscii()
    iohandler.inventory.filename = "data/matinitialize.spatialdb"
    iohandler._configure()
    db = SimpleDB()
    db.inventory.label = "material properties"
    db.inventory.iohandler = iohandler
    db._configure()

    from pylith.materials.ElasticPlaneStrain import ElasticPlaneStrain
    material = ElasticPlaneStrain()
    material.inventory.quadrature = quadrature
    material.inventory.dbProperties = db
    material.inventory.label = "my material"
    material.inventory.id = 54
    material._configure()

    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/twoelems.mesh"
    importer.inventory.coordsys = cs
    importer._configure()
    mesh = importer.read(debug=False, interpolate=False)
    
    material.preinitialize(mesh)

    # No test of result.
    return

  def _preinitialize(self):
    """
    Setup mesh and integrator and preinitialize integrator.
    """
    from spatialdata.units.Nondimensional import Nondimensional
    normalizer = Nondimensional()
    normalizer._configure()

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

    # Setup material
    from pylith.feassemble.FIATSimplex import FIATSimplex
    cell = FIATSimplex()
    cell.inventory.dimension = 2
    cell.inventory.degree = 1
    cell.inventory.order = 1
    cell._configure()
    from pylith.feassemble.Quadrature import Quadrature
    quadrature = Quadrature()
    quadrature.inventory.cell = cell
    quadrature._configure()
    
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
    iohandler = SimpleIOAscii()
    iohandler.inventory.filename = "data/elasticplanestrain.spatialdb"
    iohandler._configure()
    db = SimpleDB()
    db.inventory.label = "elastic plane strain"
    db.inventory.iohandler = iohandler
    db._configure()

    from pylith.materials.ElasticPlaneStrain import ElasticPlaneStrain
    material = ElasticPlaneStrain()
    material.inventory.label = "elastic plane strain"
    material.inventory.id = 0
    material.inventory.dbProperties = db
    material.inventory.quadrature = quadrature
    material._configure()
    
    from pylith.meshio.OutputMatElastic import OutputMatElastic
    material.output = OutputMatElastic()
    material.output._configure()
    material.output.writer._configure()

    # Setup integrator
    integrator = ElasticityImplicitLgDeform()
    integrator.preinitialize(mesh, material)
    return (mesh, integrator)

  def test(self):
    """
    Test GenSimpleDBApp with 1-D data in 2-D space.
    """
    from spatialdata.spatialdb.generator.GenSimpleDBApp import GenSimpleDBApp
    from spatialdata.geocoords.CSCart import CSCart
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    
    app = GenSimpleDBApp()
    app.run()

    # Test write using query
    db = SimpleDB()
    db.inventory.iohandler.inventory.filename = \
        "data/gen1Din2D.spatialdb"
    db.inventory.iohandler._configure()
    db.inventory.label = "test"
    db.inventory.queryType = "nearest"
    db._configure()

    qlocs = numpy.array( [[-2.0,  2.0],
                          [ 3.0, -4.0],
                          [ 0.0,  0.7]],
                         numpy.float64)
    
    dataE = numpy.array( [[-0.15, 3.45],
                          [2.4, 6.4],
                          [-0.6, 3.45]], numpy.float64)
    errE = [0, 0, 0]

    from spatialdata.geocoords.CSCart import CSCart
    cs = CSCart()
    cs.inventory.spaceDim = 2
    cs._configure()

    db.open()
    db.queryVals(["two", "one"])
    data = numpy.zeros(dataE.shape, dtype=numpy.float64)
    err = []
    nlocs = qlocs.shape[0]
    for i in xrange(nlocs):
      e = db.query(data[i,:], qlocs[i,:], cs)
      err.append(e)
    db.close()    

    self.assertEqual(len(errE), len(err))
    for vE, v in zip(errE, err):
      self.assertEqual(vE, v)

    self.assertEqual(len(dataE.shape), len(data.shape))
    for dE, d in zip(dataE.shape, data.shape):
      self.assertEqual(dE, d)
    for vE, v in zip(numpy.reshape(dataE, -1), numpy.reshape(data, -1)):
      self.assertAlmostEqual(vE, v, 6)

    return

 def setUp(self):
   from spatialdata.spatialdb.SimpleDB import SimpleDB
   db = SimpleDB()
   db.inventory.label = "test"
   db.inventory.queryType = "nearest"
   db.inventory.iohandler.inventory.filename = "data/simple.spatialdb"
   db.inventory.iohandler._configure()
   db._configure()
   self._db = db
   return

  def testDBInitialState(self):
    """
    Test dbInitialState().
    """
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
    iohandler = SimpleIOAscii()
    iohandler.inventory.filename = "data/matinitialize.spatialdb"
    iohandler._configure()
    db = SimpleDB()
    db.inventory.label = "material properties"
    db.inventory.iohandler = iohandler
    db._configure()

    self.material.dbInitialState(db)

    # No test of result.
    return

  def testDBProperties(self):
    """
    Test dbProperties().
    """
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
    iohandler = SimpleIOAscii()
    iohandler.inventory.filename = "data/staticfriction.spatialdb"
    iohandler._configure()
    db = SimpleDB()
    db.inventory.label = "friction properties"
    db.inventory.iohandler = iohandler
    db._configure()

    self.friction.dbProperties(db)

    # No test of result.
    return

  def _initialize(self):
    """
    Initialize DirichletBC boundary condition.
    """
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    dbInitial = SimpleDB()
    dbInitial.inventory.label = "TestDirichletBC tri3"
    dbInitial.inventory.iohandler.inventory.filename = "data/tri3_disp.spatialdb"
    dbInitial.inventory.iohandler._configure()
    dbInitial._configure()

    from pylith.bc.DirichletBC import DirichletBC
    bc = DirichletBC()
    bc.inventory.label = "bc"
    bc.inventory.bcDOF = [1]
    bc.inventory.dbInitial = dbInitial
    bc._configure()

    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()
    mesh = importer.read(debug=False, interpolate=False)
    
    bc.preinitialize(mesh)
    bc.initialize(totalTime=0.0, numTimeSteps=1, normalizer=normalizer)

    # Setup field
    from pylith.topology.Field import Field
    field = Field(mesh)
    field.newSection(field.VERTICES_FIELD, cs.spaceDim())
    
    return (mesh, bc, field)

  def testDBInitialStrain(self):
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
    iohandler = SimpleIOAscii()
    iohandler.inventory.filename = "data/matinitialize.spatialdb"
    iohandler._configure()
    db = SimpleDB()
    db.inventory.label = "material properties"
    db.inventory.iohandler = iohandler
    db._configure()

    material = self.material
    material.inventory.label = "Elastic material"
    material.inventory.dbInitialStrain = db
    material.inventory.useInitialStrain = True
    material._configure()

    # No test of result.
    return

  def test_initialize(self):
    """
    Test initialize().
    """
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii

    ioFinalSlip = SimpleIOAscii()
    ioFinalSlip.inventory.filename = "finalslip.spatialdb"
    ioFinalSlip._configure()
    dbFinalSlip = SimpleDB()
    dbFinalSlip.inventory.iohandler = ioFinalSlip
    dbFinalSlip.inventory.label = "final slip"
    dbFinalSlip._configure()
    
    ioSlipTime = SimpleIOAscii()
    ioSlipTime.inventory.filename = "sliptime.spatialdb"
    ioSlipTime._configure()
    dbSlipTime = SimpleDB()
    dbSlipTime.inventory.iohandler = ioSlipTime
    dbSlipTime.inventory.label = "slip time"
    dbSlipTime._configure()
    
    slipFn = StepSlipFn()
    slipFn.inventory.dbSlip = dbFinalSlip
    slipFn.inventory.dbSlipTime = dbSlipTime
    slipFn._configure()
    slipFn.preinitialize()
    slipFn.verifyConfiguration()
    slipFn.initialize()
    return

  def test_configure(self):
    """
    Test initialize().
    """
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
    from pyre.units.time import second

    ioInitial = SimpleIOAscii()
    ioInitial.inventory.filename = "tri3_initialtractions.spatialdb"
    ioInitial._configure()
    dbInitial = SimpleDB()
    dbInitial.inventory.iohandler = ioInitial
    dbInitial.inventory.label = "initial tractions"
    dbInitial._configure()
    
    ioChange = SimpleIOAscii()
    ioChange.inventory.filename = "tri3_changetractions.spatialdb"
    ioChange._configure()
    dbChange = SimpleDB()
    dbChange.inventory.iohandler = ioChange
    dbChange.inventory.label = "traction change"
    dbChange._configure()
    
    tract = TractPerturbation()
    tract.inventory.dbInitial = dbInitial
    tract.inventory.dbChange = dbChange
    tract._configure()
    return

  def test_configure(self):
    """
    Test constructor.
    """
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    db = SimpleDB()
    db.inventory.label = "simple database"
    db._configure()

    from spatialdata.spatialdb.TimeHistory import TimeHistory
    th = TimeHistory()
    th._configure()

    from pylith.bc.DirichletBC import DirichletBC
    bc = DirichletBC()
    bc.inventory.label = "abc"
    bc.inventory.dbInitial = db
    bc.inventory.dbRate = db
    bc.inventory.dbChange = db
    bc.inventory.thChange = th    
    bc._configure()
    return

  def test_initialize(self):
    """
    Test initialize().
    """
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
    from pyre.units.time import second

    ioFinalSlip = SimpleIOAscii()
    ioFinalSlip.inventory.filename = "finalslip.spatialdb"
    ioFinalSlip._configure()
    dbFinalSlip = SimpleDB()
    dbFinalSlip.inventory.iohandler = ioFinalSlip
    dbFinalSlip.inventory.label = "final slip"
    dbFinalSlip._configure()
    
    ioSlipTime = SimpleIOAscii()
    ioSlipTime.inventory.filename = "sliptime.spatialdb"
    ioSlipTime._configure()
    dbSlipTime = SimpleDB()
    dbSlipTime.inventory.iohandler = ioSlipTime
    dbSlipTime.inventory.label = "slip time"
    dbSlipTime._configure()
    
    from pylith.faults.StepSlipFn import StepSlipFn
    slipfn = StepSlipFn()
    slipfn.inventory.dbSlip = dbFinalSlip
    slipfn.inventory.dbSlipTime = dbSlipTime
    slipfn._configure()

    eqsrc = EqKinSrc()
    eqsrc.inventory.originTime = 5.3*second
    eqsrc.inventory.slipfn = slipfn
    eqsrc._configure()
    eqsrc.preinitialize()
    eqsrc.verifyConfiguration()
    eqsrc.initialize()
    return

  def test_initialize(self):
    """
    Test initialize().
    """
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
    from spatialdata.spatialdb.TimeHistory import TimeHistory

    ioFinalSlip = SimpleIOAscii()
    ioFinalSlip.inventory.filename = "finalslip.spatialdb"
    ioFinalSlip._configure()
    dbFinalSlip = SimpleDB()
    dbFinalSlip.inventory.iohandler = ioFinalSlip
    dbFinalSlip.inventory.label = "slip amplitude"
    dbFinalSlip._configure()
    
    ioSlipTime = SimpleIOAscii()
    ioSlipTime.inventory.filename = "sliptime.spatialdb"
    ioSlipTime._configure()
    dbSlipTime = SimpleDB()
    dbSlipTime.inventory.iohandler = ioSlipTime
    dbSlipTime.inventory.label = "slip time"
    dbSlipTime._configure()
    
    dbTimeHistory = TimeHistory()
    dbTimeHistory.inventory.filename = "slipfn.timedb"
    dbTimeHistory.inventory.label = "time history"
    dbTimeHistory._configure()
    
    slipFn = TimeHistorySlipFn()
    slipFn.inventory.dbslip = dbFinalSlip
    slipFn.inventory.dbSlipTime = dbSlipTime
    slipFn.inventory.dbTimeHistory = dbTimeHistory
    slipFn._configure()
    slipFn.preinitialize()
    slipFn.verifyConfiguration()
    slipFn.initialize()
    return

  def test_write(self):
    """
    Test write().
    """
    # Database info
    cs = CSCart()
    cs.initialize()

    filename = "data/test.spatialdb"
    data = {'points': numpy.array( [ [1.0, 2.0, 3.0],
                                   [0.5, 3.0, -3.0]], numpy.float64),
            'coordsys': cs,
            'data_dim': 1,
            'values': [{'name': "One",
                        'units': "m",
                        'data': numpy.array( [2.0, 8.0], numpy.float64)},
                       {'name': "Two",
                        'units': "m",
                        'data': numpy.array( [-2.0, 3.0], numpy.float64)}]}            
    dataDim = 1
            
    qlocs = numpy.array( [[0.875, 2.25, 1.5],
                          [0.6, 2.8, -1.8],
                          [1.0, 2.0, 3.0]],
                         numpy.float64)    
    valsE = numpy.array( [[-0.75, 3.5],
                          [2.0, 6.8],
                          [-2.0, 2.0]], numpy.float64)
    errE = [0, 0, 0]
  
  
    # Write database
    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
    writer = SimpleIOAscii()
    writer.inventory.filename = filename
    writer._configure()
    writer.write(data)

    # Test write using query
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    db = SimpleDB()
    db.inventory.label = "test"
    db.inventory.queryType = "linear"
    db.inventory.iohandler.inventory.filename = filename
    db.inventory.iohandler._configure()
    db._configure()

    db.open()
    db.queryVals(["two", "one"])
    vals = numpy.zeros(valsE.shape, dtype=numpy.float64)
    err = []
    nlocs = qlocs.shape[0]
    for i in xrange(nlocs):
      e = db.query(vals[i,:], qlocs[i,:], cs)
      err.append(e)
    db.close()    
    
    self.assertEqual(len(valsE.shape), len(vals.shape))
    for dE, d in zip(valsE.shape, vals.shape):
      self.assertEqual(dE, d)
    for vE, v in zip(numpy.reshape(valsE, -1), numpy.reshape(vals, -1)):
      self.assertAlmostEqual(vE, v, 6)

    return

  def _initialize(self):
    """
    Initialize AbsorbingDampers boundary condition.
    """
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    db = SimpleDB()
    db.inventory.label = "TestAbsorbingDampers tri3"
    db.inventory.iohandler.inventory.filename = \
        "data/elasticplanestrain.spatialdb"
    db.inventory.iohandler._configure()
    db._configure()

    from pylith.feassemble.FIATSimplex import FIATSimplex
    cell = FIATSimplex()
    cell.inventory.dimension = 1
    cell.inventory.degree = 1
    cell.inventory.order = 1
    cell._configure()
    from pylith.feassemble.Quadrature import Quadrature
    quadrature = Quadrature()
    quadrature.inventory.cell = cell
    quadrature._configure()

    from pylith.bc.AbsorbingDampers import AbsorbingDampers
    bc = AbsorbingDampers()
    bc.inventory.quadrature = quadrature
    bc.inventory.db = db
    bc.inventory.id = 0
    bc.inventory.label = "bc"
    bc._configure()

    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()
    mesh = importer.read(debug=False, interpolate=False)
    
    bc.preinitialize(mesh)
    bc.initialize(totalTime=0.0, numTimeSteps=1, normalizer=normalizer)
    bc.timeStep(0.01)

    # Setup fields
    from pylith.topology.SolutionFields import SolutionFields
    fields = SolutionFields(mesh)
    fields.add("residual", "residual")
    fields.add("dispIncr(t->t+dt)", "displacement")
    fields.add("disp(t)", "displacement")
    fields.add("disp(t-dt)", "displacement")
    fields.add("velocity(t)", "velocity")
    fields.solutionName("dispIncr(t->t+dt)")

    residual = fields.get("residual")
    residual.newSection(residual.VERTICES_FIELD, cs.spaceDim())
    residual.allocate()
    residual.zero()

    fields.copyLayout("residual")
    
    return (mesh, bc, fields)

  def _initialize(self):
    """
    Initialize fault.
    """
    dt = 2.4
    
    from spatialdata.units.Nondimensional import Nondimensional
    normalizer = Nondimensional()
    normalizer._configure()

    # Setup mesh
    cs = CSCart()
    cs.inventory.spaceDim = 2
    cs._configure()

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

    # Setup quadrature
    from pylith.feassemble.FIATSimplex import FIATSimplex
    cell = FIATSimplex()
    cell.inventory.dimension = 1
    cell.inventory.degree = 1
    cell.inventory.order = 1
    cell._configure()
    from pylith.feassemble.Quadrature import Quadrature
    quadrature = Quadrature()
    quadrature.inventory.cell = cell
    quadrature._configure()

    # Setup rupture info
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
    ioTractions = SimpleIOAscii()
    ioTractions.inventory.filename = "data/tri3_initialtractions.spatialdb"
    ioTractions._configure()
    dbTractions = SimpleDB()
    dbTractions.inventory.iohandler = ioTractions
    dbTractions.inventory.label = "initial tractions"
    dbTractions._configure()
    from pylith.faults.TractPerturbation import TractPerturbation
    tract = TractPerturbation()
    tract.inventory.dbInitial = dbTractions
    tract._configure()

    ioFriction = SimpleIOAscii()
    ioFriction.inventory.filename = "data/tri3_staticfriction.spatialdb"
    ioFriction._configure()
    dbFriction = SimpleDB()
    dbFriction.inventory.iohandler = ioFriction
    dbFriction.inventory.label = "friction"
    dbFriction._configure()
    
    from pylith.friction.StaticFriction import StaticFriction
    friction = StaticFriction()
    friction.inventory.label = "Static friction"
    friction.inventory.dbProperties = dbFriction
    friction._configure()

    # Setup fault
    fault = FaultCohesiveDyn()
    fault.inventory.output.inventory.writer._configure()
    fault.inventory.output._configure()
    fault.inventory.matId = 10
    fault.inventory.faultLabel = "fault"
    fault.inventory.upDir = [0, 0, 1]
    fault.inventory.faultQuadrature = quadrature
    fault.inventory.tract = tract
    fault.inventory.friction = friction
    fault._configure()

    nvertices = fault.numVerticesNoMesh(mesh)
    firstFaultVertex = 0
    firstLagrangeVertex = nvertices
    firstFaultCell      = 2*nvertices
    fault.adjustTopology(mesh, firstFaultVertex, firstLagrangeVertex,
                         firstFaultCell)
    from pylith.topology.topology import MeshOps_nondimensionalize
    MeshOps_nondimensionalize(mesh, normalizer)

    fault.preinitialize(mesh)
    fault.timeStep(dt)
    fault.verifyConfiguration()
    from pyre.units.time import s
    fault.initialize(totalTime=0.0*s, numTimeSteps=1, normalizer=normalizer)

    # Setup fields
    from pylith.topology.SolutionFields import SolutionFields
    fields = SolutionFields(mesh)
    fields.add("residual", "residual")
    fields.add("dispIncr(t->t+dt)", "displacement_increment")
    fields.add("disp(t)", "displacement")
    fields.add("velocity(t)", "velocity")
    fields.solutionName("dispIncr(t->t+dt)")

    residual = fields.get("residual")
#.........这里部分代码省略.........

  def _initialize(self):
    """
    Initialize fault.
    """
    dt = 2.4
    
    from spatialdata.units.Nondimensional import Nondimensional
    normalizer = Nondimensional()
    normalizer._configure()

    # Setup mesh
    cs = CSCart()
    cs.inventory.spaceDim = 2
    cs._configure()

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

    # Setup quadrature
    from pylith.feassemble.FIATSimplex import FIATSimplex
    cell = FIATSimplex()
    cell.inventory.dimension = 1
    cell.inventory.degree = 1
    cell.inventory.order = 1
    cell._configure()
    from pylith.feassemble.Quadrature import Quadrature
    quadrature = Quadrature()
    quadrature.inventory.cell = cell
    quadrature._configure()

    # Setup earthquake source
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
    ioFinalSlip = SimpleIOAscii()
    ioFinalSlip.inventory.filename = "data/tri3_finalslip.spatialdb"
    ioFinalSlip._configure()
    dbFinalSlip = SimpleDB()
    dbFinalSlip.inventory.iohandler = ioFinalSlip
    dbFinalSlip.inventory.label = "final slip"
    dbFinalSlip._configure()
    
    ioSlipTime = SimpleIOAscii()
    ioSlipTime.inventory.filename = "data/tri3_sliptime.spatialdb"
    ioSlipTime._configure()
    dbSlipTime = SimpleDB()
    dbSlipTime.inventory.iohandler = ioSlipTime
    dbSlipTime.inventory.label = "slip time"
    dbSlipTime._configure()
    
    from pylith.faults.StepSlipFn import StepSlipFn
    slipfn = StepSlipFn()
    slipfn.inventory.dbSlip = dbFinalSlip
    slipfn.inventory.dbSlipTime = dbSlipTime
    slipfn._configure()

    # Setup fault
    fault = FaultCohesiveKin()
    fault.inventory.output.inventory.writer._configure()
    fault.inventory.output._configure()
    fault.inventory.matId = 10
    fault.inventory.faultLabel = "fault"
    fault.inventory.upDir = [0, 0, 1]
    fault.inventory.faultQuadrature = quadrature
    fault._configure()
    eqsrc = fault.eqsrcs.components()[0]
    eqsrc.inventory.originTime = 1.23*second
    eqsrc.inventory.slipfn = slipfn
    eqsrc._configure()

    nvertices = fault.numVerticesNoMesh(mesh)
    firstFaultVertex = 0
    firstLagrangeVertex = nvertices
    firstFaultCell      = 2*nvertices
    fault.adjustTopology(mesh, firstFaultVertex, firstLagrangeVertex,
                         firstFaultCell)
    fault.preinitialize(mesh)
    fault.timeStep(dt)
    fault.verifyConfiguration()
    from pyre.units.time import s
    fault.initialize(totalTime=0.0*s, numTimeSteps=1, normalizer=normalizer)

    # Setup fields
    from pylith.topology.SolutionFields import SolutionFields
    fields = SolutionFields(mesh)
    fields.add("residual", "residual")
    fields.add("dispIncr(t->t+dt)", "displacement_increment")
    fields.add("disp(t)", "displacement")
    fields.solutionName("dispIncr(t->t+dt)")

    residual = fields.get("residual")
    residual.subfieldAdd("displacement", cs.spaceDim(), residual.VECTOR)
    residual.subfieldAdd("lagrange_multiplier", cs.spaceDim(), residual.VECTOR)
    residual.subfieldsSetup()
    residual.setupSolnChart()
    residual.setupSolnDof(cs.spaceDim())
    fault.setupSolnDof(residual)
#.........这里部分代码省略.........

  def _initialize(self):
    """
    Initialize fault.
    """
    dt = 2.4
    
    from spatialdata.units.Nondimensional import Nondimensional
    normalizer = Nondimensional()
    normalizer._configure()

    # Setup mesh
    cs = CSCart()
    cs.inventory.spaceDim = 2
    cs._configure()

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

    # Setup quadrature
    from pylith.feassemble.FIATSimplex import FIATSimplex
    cell = FIATSimplex()
    cell.inventory.dimension = 1
    cell.inventory.degree = 1
    cell.inventory.order = 1
    cell._configure()
    from pylith.feassemble.Quadrature import Quadrature
    quadrature = Quadrature()
    quadrature.inventory.cell = cell
    quadrature._configure()

    # Setup impulses
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
    ioImpulseAmp = SimpleIOAscii()
    ioImpulseAmp.inventory.filename = "data/tri3_impulses.spatialdb"
    ioImpulseAmp._configure()
    dbImpulseAmp = SimpleDB()
    dbImpulseAmp.inventory.iohandler = ioImpulseAmp
    dbImpulseAmp.inventory.label = "impulse amplitude"
    dbImpulseAmp._configure()
    
    # Setup fault
    fault = FaultCohesiveImpulses()
    fault.inventory.output.inventory.writer._configure()
    fault.inventory.output._configure()
    fault.inventory.matId = 10
    fault.inventory.faultLabel = "fault"
    fault.inventory.upDir = [0, 0, 1]
    fault.inventory.faultQuadrature = quadrature
    fault.inventory.dbImpulseAmp = dbImpulseAmp
    fault._configure()

    nvertices = fault.numVerticesNoMesh(mesh)
    firstFaultVertex = 0
    firstLagrangeVertex = nvertices
    firstFaultCell      = 2*nvertices
    fault.adjustTopology(mesh, firstFaultVertex, firstLagrangeVertex,
                         firstFaultCell)
    fault.preinitialize(mesh)
    fault.timeStep(dt)
    fault.verifyConfiguration()
    from pyre.units.time import s
    fault.initialize(totalTime=0.0*s, numTimeSteps=1, normalizer=normalizer)

    # Setup fields
    from pylith.topology.SolutionFields import SolutionFields
    fields = SolutionFields(mesh)
    fields.add("residual", "residual")
    fields.add("dispIncr(t->t+dt)", "displacement_increment")
    fields.add("disp(t)", "displacement")
    fields.solutionName("dispIncr(t->t+dt)")
    residual = fields.get("residual")
    residual.newSection(residual.VERTICES_FIELD, cs.spaceDim())
    residual.allocate()
    residual.zero()
    fields.copyLayout("residual")
    
    return (mesh, fault, fields)

def getCellSizeDB(points):

    from spatialdata.geocoords.CSCart import CSCart
    from spatialdata.spatialdb.SimpleDB import SimpleDB
    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii

    # Coordinate system for mesh (must match coordsys in ExodusII file)
    cs = CSCart()
    cs._configure()
    cs.initialize()

    # Spatial database with physical properties (Vs)
    dbIO = SimpleIOAscii()
    dbIO.inventory.filename = filenameDB
    dbIO._configure()
    db = SimpleDB()
    db.inventory.iohandler = dbIO
    db.inventory.label = "Physical properties"
    db.inventory.queryType = "linear"
    db._configure()
    
    (npoints, spacedim) = points.shape

    # Query database
    db.open()
    db.queryVals(["vs"])
    data = numpy.zeros((npoints, 1), dtype=numpy.float64)
    err = numpy.zeros((npoints,), dtype=numpy.int32)
    db.multiquery(data, err, points, cs)
    db.close()

    vs = data[:,0]
    cellSize = minPeriod*vs / 10.0
    return cellSize