Python SolutionFields.SolutionFields类代码示例

  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 setUp(self):
    from pylith.meshio.MeshIOAscii import MeshIOAscii
    iohandler = MeshIOAscii()
    filename = "data/twohex8.txt"
    
    from spatialdata.units.Nondimensional import Nondimensional
    normalizer = Nondimensional()
    normalizer._configure()

    from spatialdata.geocoords.CSCart import CSCart
    iohandler.inventory.filename = filename
    iohandler.inventory.coordsys = CSCart()
    iohandler._configure()
    mesh = iohandler.read(debug=False, interpolate=False)

    from pylith.topology.SolutionFields import SolutionFields
    fields = SolutionFields(mesh)

    name = "disp(t)"
    fields.add(name, "displacement")
    fields.solutionName(name)
    field = fields.get(name)
    field.subfieldAdd("displacement", mesh.dimension(), field.VECTOR)
    field.subfieldsSetup()
    field.newSection(field.VERTICES_FIELD, mesh.dimension())
    field.allocate()

    self.mesh = mesh
    self.fields = fields
    self.normalizer = normalizer
    return

  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)
    
    self.fields = SolutionFields(self.mesh)
    return

  def _initialize(self, mesh, integrator):
    """
    Initialize integrator.
    """
    dt = 2.3
    
    from spatialdata.units.Nondimensional import Nondimensional
    normalizer = Nondimensional()
    normalizer._configure()

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

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

    residual = fields.get("residual")
    spaceDim = mesh.coordsys().spaceDim()
    lengthScale = normalizer.lengthScale()
    residual.subfieldAdd("displacement", spaceDim, residual.VECTOR, lengthScale.value);
    residual.subfieldAdd("lagrange_multiplier", spaceDim, residual.VECTOR);
    
    residual.subfieldsSetup();
    residual.setupSolnChart();
    residual.setupSolnDof(spaceDim);
    residual.allocate();
    residual.zeroAll();
    fields.copyLayout("residual")

    return fields

  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)

class TestSolutionFields(unittest.TestCase):
  """
  Unit testing of SolutionFields 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)
    
    self.fields = SolutionFields(self.mesh)
    return


  def test_constructor(self):
    """
    Test constructor.
    """
    return


  def test_solutionName(self):
    """
    Test mesh().
    """
    fields = self.fields
    fields.add("field A", "A");
    fields.add("field B", "B");
    fields.add("field C", "C");
    
    fields.solutionName("field B")
    return


  def test_solution(self):
    """
    Test solution().
    """
    fields = self.fields
    fields.add("field A", "A");
    fields.add("field B", "B");
    fields.add("field C", "C");
    
    fields.solutionName("field B")
    solution = self.fields.solution()
    return


  def test_fieldAdd(self):
    """
    Test fieldAdd().
    """
    fields = self.fields
    fields.add("field A", "A");
    fields.add("field B", "B");

    helper_fieldAdd(fields)
    fieldA = fields.get("field A")
    fieldB = fields.get("field B")
    fieldA.allocate()
    fieldB.allocate()
    fieldA.copy(fieldB)
    return

  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, mesh, integrator):
    """
    Initialize integrator.
    """
    dt = 2.3
    
    from spatialdata.units.Nondimensional import Nondimensional
    normalizer = Nondimensional()
    normalizer._configure()

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

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

    residual = fields.get("residual")
    residual.newSection(residual.VERTICES_FIELD, mesh.coordsys().spaceDim())
    residual.allocate()
    fields.copyLayout("residual")

    residual.zero()
    return 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 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 _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

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