[cig-commits] r7016 - in short/3D/PyLith/trunk: . pylith/feassemble pylith/problems unittests/libtests/feassemble unittests/pytests/feassemble unittests/pytests/feassemble/data

brad at geodynamics.org brad at geodynamics.org
Thu May 31 12:06:11 PDT 2007


Author: brad
Date: 2007-05-31 12:06:10 -0700 (Thu, 31 May 2007)
New Revision: 7016

Added:
   short/3D/PyLith/trunk/unittests/pytests/feassemble/data/
   short/3D/PyLith/trunk/unittests/pytests/feassemble/data/Makefile.am
   short/3D/PyLith/trunk/unittests/pytests/feassemble/data/elasticplanestrain.spatialdb
   short/3D/PyLith/trunk/unittests/pytests/feassemble/data/tri3.mesh
Modified:
   short/3D/PyLith/trunk/TODO
   short/3D/PyLith/trunk/configure.ac
   short/3D/PyLith/trunk/pylith/feassemble/Integrator.py
   short/3D/PyLith/trunk/pylith/problems/Explicit.py
   short/3D/PyLith/trunk/pylith/problems/Implicit.py
   short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
   short/3D/PyLith/trunk/unittests/pytests/feassemble/Makefile.am
   short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py
   short/3D/PyLith/trunk/unittests/pytests/feassemble/testfeassemble.py
Log:
Finished implementing Python unit tests for explicit elasticity integrator. Fixed some bugs exposed by the unit tests.

Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO	2007-05-31 17:52:01 UTC (rev 7015)
+++ short/3D/PyLith/trunk/TODO	2007-05-31 19:06:10 UTC (rev 7016)
@@ -2,38 +2,38 @@
 MAIN PRIORITIES (Brad)
 ======================================================================
 
-1. Unit tests for Integrator stuff.
+topology
+  add section->add(fieldA, fieldB)
 
-  a. Add in use of useElasticBehavior()
+ElasticityExplicit
+  use new mesh->restrict() [preallocated array]
 
-  b. Add unit test for IntegratorElasticity::calcTotalStrain
+ElasticityImplicit
+  use new mesh->restrict() [preallocated array]
 
-  c. ElasticityExplicit
-    i. C++ unit tests
-      Need help from Matt to fix PETSc matrix update problem.
-    ii. Python unit tests
-      (1) setMesh()
-      (2) initQuadrature()
-      (3) timeStep()
-      (4) stableTimeStep()
-      (5) integrateResidual() [single test is okay, inteface is same]
-      (6) integrateJacobian() [single test is okay, inteface is same]
-      (7) updateState()
+Implicit.py
+  use section->add(fieldA, fieldB)
 
-  d. ElasticityImplicit
+Add useElasticBehavior() switches to Implicit::solveElastic().
+
+1. Unit tests for Integrator stuff.
+
+  a. Add unit test for IntegratorElasticity::calcTotalStrain
+
+  c. ElasticityImplicit
     i. C++ unit tests
       (1) timeStep()
       (2) stableTimeStep()
-      (3) integrateResidual() [Charles]
-      (4) integrateJacobian() [Charles]
+      (3) integrateResidual()
+      (4) integrateJacobian()
       (5) updateState()
     ii. Python unit tests
       (1) setMesh()
       (2) initQuadrature()
       (3) timeStep()
       (4) stableTimeStep()
-      (5) integrateResidual() [Charles]
-      (6) integrateJacobian() [Charles]
+      (5) integrateResidual()
+      (6) integrateJacobian()
       (7) updateState()
 
 2. Add dualBasis to Quadrature.

Modified: short/3D/PyLith/trunk/configure.ac
===================================================================
--- short/3D/PyLith/trunk/configure.ac	2007-05-31 17:52:01 UTC (rev 7015)
+++ short/3D/PyLith/trunk/configure.ac	2007-05-31 19:06:10 UTC (rev 7016)
@@ -238,6 +238,7 @@
 		unittests/pytests/faults/Makefile
 		unittests/pytests/faults/data/Makefile
 		unittests/pytests/feassemble/Makefile
+		unittests/pytests/feassemble/data/Makefile
 		unittests/pytests/materials/Makefile
 		unittests/pytests/materials/data/Makefile
 		unittests/pytests/meshio/Makefile

Modified: short/3D/PyLith/trunk/pylith/feassemble/Integrator.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/Integrator.py	2007-05-31 17:52:01 UTC (rev 7015)
+++ short/3D/PyLith/trunk/pylith/feassemble/Integrator.py	2007-05-31 19:06:10 UTC (rev 7016)
@@ -94,7 +94,8 @@
     Integrate contributions to residual term at time t.
     """
     assert(None != self.cppHandle)
-    self.cppHandle.integrateResidual(residual, fields, self.mesh.cppHandle)
+    self.cppHandle.integrateResidual(residual, fields.cppHandle,
+                                     self.mesh.cppHandle)
     return
 
 
@@ -112,7 +113,8 @@
     Integrate contributions to Jacobian term at time t.
     """
     assert(None != self.cppHandle)
-    self.cppHandle.integrateJacobian(jacobian, fields, self.mesh.cppHandle)
+    self.cppHandle.integrateJacobian(jacobian, fields.cppHandle,
+                                     self.mesh.cppHandle)
     return
 
 

Modified: short/3D/PyLith/trunk/pylith/problems/Explicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Explicit.py	2007-05-31 17:52:01 UTC (rev 7015)
+++ short/3D/PyLith/trunk/pylith/problems/Explicit.py	2007-05-31 19:06:10 UTC (rev 7016)
@@ -84,7 +84,7 @@
     petsc.mat_setzero(self.jacobian)
     for integrator in self.integrators:
       integrator.timeStep(dt)
-      integrator.integrateJacobian(self.jacobian, self.fields.cppHandle)
+      integrator.integrateJacobian(self.jacobian, self.fields)
     petsc.mat_assemble(self.jacobian)
 
     self.solver.initialize(mesh, self.fields.getReal("dispTpdt"))
@@ -119,7 +119,7 @@
       petsc.mat_setzero(self.jacobian)
       for integrator in self.integrators:
         integrator.timeStep(dt)
-        integrator.integrateJacobian(self.jacobian, self.fields.cppHandle)
+        integrator.integrateJacobian(self.jacobian, self.fields)
       petsc.mat_assemble(self.jacobian)
     return
 
@@ -134,7 +134,7 @@
     bindings.zeroRealSection(residual)
     for integrator in self.integrators:
       integrator.timeStep(dt)
-      integrator.integrateResidual(residual, self.fields.cppHandle)
+      integrator.integrateResidual(residual, self.fields)
 
     self._info.log("Solving equations.")
     self.solver.solve(self.fields.getReal("dispTpdt"), self.jacobian, residual)

Modified: short/3D/PyLith/trunk/pylith/problems/Implicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Implicit.py	2007-05-31 17:52:01 UTC (rev 7015)
+++ short/3D/PyLith/trunk/pylith/problems/Implicit.py	2007-05-31 19:06:10 UTC (rev 7016)
@@ -139,7 +139,7 @@
       petsc.mat_setzero(self.jacobian)
       for integrator in self.integrators:
         integrator.timeStep(dt)
-        integrator.integrateJacobian(self.jacobian, self.fields.cppHandle)
+        integrator.integrateJacobian(self.jacobian, self.fields)
       petsc.mat_assemble(self.jacobian)
     # Put in loop over integrators to see if stiffness needs
     # reforming.  If so, then reform it.
@@ -156,7 +156,7 @@
     bindings.zeroRealSection(residual)
     for integrator in self.integrators:
       integrator.timeStep(dt)
-      integrator.integrateResidual(residual, self.fields.cppHandle)
+      integrator.integrateResidual(residual, self.fields)
 
     self._info.log("Solving equations.")
     self.solver.solve(self.fields.getReal("dispIncr"), self.jacobian, residual)
@@ -217,9 +217,9 @@
     petsc.mat_setzero(self.jacobian)
     for integrator in self.integrators:
       integrator.timeStep(dt)
-      integrator.integrateJacobian(self.jacobian, self.fields.cppHandle)
+      integrator.integrateJacobian(self.jacobian, self.fields)
       integrator.integrateResidual(self.fields.getReal("dispT"),
-                                   self.fields.cppHandle)
+                                   self.fields)
     import pylith.utils.petsc as petsc
     petsc.mat_assemble(self.jacobian)
 

Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc	2007-05-31 17:52:01 UTC (rev 7015)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityExplicit.cc	2007-05-31 19:06:10 UTC (rev 7016)
@@ -144,7 +144,6 @@
   ElasticityExplicit integrator;
   topology::FieldsManager fields(mesh);
   _initialize(&mesh, &integrator, &fields);
-  mesh->getFactory()->clear();
 
   const ALE::Obj<pylith::real_section_type>& dispTpdt = 
     fields.getReal("dispTpdt");
@@ -235,6 +234,7 @@
       c_iter != cells->end();
       ++c_iter)
     (*mesh)->setValue(labelMaterials, *c_iter, _data->matId);
+  (*mesh)->getFactory()->clear(); // clear numberings
 
   // Setup quadrature
   _quadrature->initialize(_data->basis, _data->basisDeriv, _data->quadPts,

Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/Makefile.am	2007-05-31 17:52:01 UTC (rev 7015)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/Makefile.am	2007-05-31 19:06:10 UTC (rev 7016)
@@ -13,6 +13,8 @@
 subpackage = feassemble
 include $(top_srcdir)/subpackage.am
 
+SUBDIRS = data
+
 TESTS = testfeassemble.py
 
 check_SCRIPTS = testfeassemble.py

Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py	2007-05-31 17:52:01 UTC (rev 7015)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityExplicit.py	2007-05-31 19:06:10 UTC (rev 7016)
@@ -17,6 +17,8 @@
 import unittest
 from pylith.feassemble.ElasticityExplicit import ElasticityExplicit
 
+from spatialdata.geocoords.CSCart import CSCart
+
 # ----------------------------------------------------------------------
 class TestElasticityExplicit(unittest.TestCase):
   """
@@ -56,7 +58,18 @@
     """
     Test setMesh().
     """
-    raise NotImplementedError("Unit test not implemented.")
+    cs = CSCart()
+    cs.spaceDim = 2
+    
+    from pylith.meshio.MeshIOAscii import MeshIOAscii
+    importer = MeshIOAscii()
+    importer.filename = "data/tri3.mesh"
+    importer.coordsys = cs
+    mesh = importer.read(debug=False, interpolate=False)
+
+    integrator = ElasticityExplicit()
+    integrator.setMesh(mesh)
+    self.assertEqual(mesh, integrator.mesh)
     return
 
   
@@ -84,37 +97,134 @@
     return
 
   
-  def test_integrateResidual(self):
+  def test_needNewJacobian(self):
     """
-    Test integrateResidual().
+    Test needNewJacobian().
     """
-    raise NotImplementedError("Unit test not implemented.")
+    integrator = ElasticityExplicit()
+    self.assertEqual(False, integrator.needNewJacobian())
     return
 
   
-  def test_needNewJacobian(self):
+  def test_integrateResidual(self):
     """
-    Test needNewJacobian().
+    Test integrateResidual().
+
+    WARNING: This is not a rigorous test of integrateResidual() because we
+    neither set the input fields or verify the results.
     """
-    integrator = ElasticityExplicit()
-    self.assertEqual(False, integrator.needNewJacobian())
+    (mesh, integrator, fields) = self._initialize()
+
+    residual = fields.getReal("residual")
+    integrator.integrateResidual(residual, fields)
+
+    # 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.
     """
-    raise NotImplementedError("Unit test not implemented.")
+    (mesh, integrator, fields) = self._initialize()
+
+    jacobian = mesh.createMatrix(fields.getReal("residual"))
+    import pylith.utils.petsc as petsc
+    petsc.mat_setzero(jacobian)
+    integrator.integrateJacobian(jacobian, fields)
+
+    # We should really add something here to check to make sure things
+    # actually initialized correctly    
     return
 
   
   def test_updateState(self):
     """
     Test updateState().
+
+    WARNING: This is not a rigorous test of updateState() because we
+    neither set the input fields or verify the results.
     """
-    raise NotImplementedError("Unit test not implemented.")
+    (mesh, integrator, fields) = self._initialize()
+
+    dispT = fields.getReal("dispT")
+    integrator.updateState(dispT)
+
+    # We should really add something here to check to make sure things
+    # actually initialized correctly    
     return
   
 
+  # PRIVATE METHODS ////////////////////////////////////////////////////
+
+  def _initialize(self):
+    """
+    Initialize integrator.
+    """
+    from pyre.units.time import second
+    dt = 1.0*second
+    
+    # Setup mesh
+    cs = CSCart()
+    cs.spaceDim = 2
+    from pylith.meshio.MeshIOAscii import MeshIOAscii
+    importer = MeshIOAscii()
+    importer.filename = "data/tri3.mesh"
+    importer.coordsys = cs
+    mesh = importer.read(debug=False, interpolate=False)
+
+    # Setup material
+    from pylith.feassemble.FIATSimplex import FIATSimplex
+    cell = FIATSimplex()
+    cell.shape = "triangle"
+    cell.degree = 1
+    cell.order = 1
+    from pylith.feassemble.quadrature.Quadrature2D import Quadrature2D
+    quadrature = Quadrature2D()
+    quadrature.cell = cell
+    
+    from spatialdata.spatialdb.SimpleDB import SimpleDB
+    from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
+    iohandler = SimpleIOAscii()
+    iohandler.filename = "data/elasticplanestrain.spatialdb"
+    db = SimpleDB()
+    db.label = "elastic plane strain"
+    db.iohandler = iohandler
+
+    from pylith.materials.ElasticPlaneStrain import ElasticPlaneStrain
+    material = ElasticPlaneStrain()
+    material.id = 0
+    material.label = "elastic plane strain"
+    material.db = db
+    material.quadrature = quadrature
+
+    # Setup integrator
+    integrator = ElasticityExplicit()
+    integrator.setMesh(mesh)
+    integrator.initQuadrature(material.quadrature)
+    integrator.initMaterial(mesh, material)
+    integrator.timeStep(dt)
+
+    # Setup fields
+    from pylith.topology.FieldsManager import FieldsManager
+    fields = FieldsManager(mesh)
+    fields.addReal("residual")
+    fields.addReal("dispTpdt")
+    fields.addReal("dispT")
+    fields.addReal("dispTmdt")
+    fields.createHistory(["dispTpdt", "dispT", "dispTmdt"])
+    fields.setFiberDimension("residual", cs.spaceDim)
+    fields.allocate("residual")
+    fields.copyLayout("residual")
+
+    import pylith.topology.topology as bindings
+    bindings.zeroRealSection(fields.getReal("residual"))
+    
+    return (mesh, integrator, fields)
+
+
 # End of file 

Added: short/3D/PyLith/trunk/unittests/pytests/feassemble/data/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/data/Makefile.am	2007-05-31 17:52:01 UTC (rev 7015)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/data/Makefile.am	2007-05-31 19:06:10 UTC (rev 7016)
@@ -0,0 +1,30 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+#                           Brad T. Aagaard
+#                        U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+noinst_DATA = \
+	tri3.mesh \
+	elasticplanestrain.spatialdb
+
+noinst_TMP =
+
+# 'export' the input files by performing a mock install
+export_datadir = $(top_builddir)/unittests/pytests/feassemble/data
+export-data: $(noinst_DATA)
+	for f in $(noinst_DATA); do $(install_sh_DATA) $(srcdir)/$$f $(export_datadir); done
+
+BUILT_SOURCES = export-data
+
+CLEANFILES = \
+	$(export_datadir)/$(noinst_DATA)
+
+
+# End of file 

Added: short/3D/PyLith/trunk/unittests/pytests/feassemble/data/elasticplanestrain.spatialdb
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/data/elasticplanestrain.spatialdb	2007-05-31 17:52:01 UTC (rev 7015)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/data/elasticplanestrain.spatialdb	2007-05-31 19:06:10 UTC (rev 7016)
@@ -0,0 +1,14 @@
+#SPATIAL.ascii 1
+SimpleDB {
+  num-values = 3
+  value-names =  density vs vp
+  value-units =  kg/m^3  m/s  m/s
+  num-locs = 1
+  data-dim = 0
+  space-dim = 2
+  cs-data = cartesian {
+    to-meters = 1.0
+    space-dim = 2
+  }
+}
+0.0  0.0   2500.0  3464.1016151377544 6000.0

Added: short/3D/PyLith/trunk/unittests/pytests/feassemble/data/tri3.mesh
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/data/tri3.mesh	2007-05-31 17:52:01 UTC (rev 7015)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/data/tri3.mesh	2007-05-31 19:06:10 UTC (rev 7016)
@@ -0,0 +1,42 @@
+mesh = {
+  dimension = 2
+  use-index-zero = true
+  vertices = {
+    dimension = 2
+    count = 4
+    coordinates = {
+             0     -1.0  0.0
+             1      0.0 -1.0
+             2      0.0  1.0
+             3      1.0  0.0
+    }
+  }
+  cells = {
+    count = 2
+    num-corners = 3
+    simplices = {
+             0       0  1  2
+             1       1  3  2
+    }
+    material-ids = {
+             0   0
+             1   0
+    }
+  }
+  group = {
+    name = bc
+    type = vertices
+    count = 2
+    indices = {
+      1  3
+    }
+  }
+  group = {
+    name = bc2
+    type = vertices
+    count = 1
+    indices = {
+      0
+    }
+  }
+}

Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/testfeassemble.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/testfeassemble.py	2007-05-31 17:52:01 UTC (rev 7015)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/testfeassemble.py	2007-05-31 19:06:10 UTC (rev 7016)
@@ -37,7 +37,11 @@
     """
     Run the application.
     """
+    from pylith.utils.PetscManager import PetscManager
+    manager = PetscManager()
+    manager.initialize()
     unittest.TextTestRunner(verbosity=2).run(self._suite())
+    manager.finalize()
     return
 
 



More information about the cig-commits mailing list