[cig-commits] r14310 - in short/3D/PyLith/branches/pylith-swig: libsrc/topology modulesrc/topology pylith/topology unittests/libtests/topology unittests/pytests/topology
brad at geodynamics.org
brad at geodynamics.org
Thu Mar 12 16:10:02 PDT 2009
Author: brad
Date: 2009-03-12 16:10:01 -0700 (Thu, 12 Mar 2009)
New Revision: 14310
Added:
short/3D/PyLith/branches/pylith-swig/libsrc/topology/Jacobian.cc
short/3D/PyLith/branches/pylith-swig/libsrc/topology/Jacobian.hh
short/3D/PyLith/branches/pylith-swig/modulesrc/topology/Jacobian.i
short/3D/PyLith/branches/pylith-swig/pylith/topology/Jacobian.py
short/3D/PyLith/branches/pylith-swig/unittests/libtests/topology/TestJacobian.cc
short/3D/PyLith/branches/pylith-swig/unittests/libtests/topology/TestJacobian.hh
short/3D/PyLith/branches/pylith-swig/unittests/pytests/topology/TestJacobian.py
Log:
Added Jacobian to wrap PETSc sparse matrix.
Added: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Jacobian.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Jacobian.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Jacobian.cc 2009-03-12 23:10:01 UTC (rev 14310)
@@ -0,0 +1,165 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "Jacobian.hh" // implementation of class methods
+
+#include "Mesh.hh" // USES Mesh
+#include "SolutionFields.hh" // USES SolutionFields
+#include "Field.hh" // USES Field
+
+// ----------------------------------------------------------------------
+// Default constructor.
+pylith::topology::Jacobian::Jacobian(const SolutionFields& fields) :
+ _fields(fields),
+ _matrix(0)
+{ // constructor
+ const ALE::Obj<Mesh::SieveMesh>& sieveMesh = fields.mesh().sieveMesh();
+ const ALE::Obj<Mesh::RealSection>& solnSection = fields.solution().section();
+
+ _matrix = new PetscMat;
+ assert(0 != _matrix);
+ PetscErrorCode err = MeshCreateMatrix(sieveMesh, solnSection,
+ MATAIJ, _matrix);
+ if (err) {
+ PetscError(__LINE__,__FUNCT__,__FILE__,__SDIR__,err,0," ");
+ throw std::runtime_error("Could not create PETSc sparse matrix "
+ "associated with system Jacobian.");
+ } // if
+
+
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor.
+pylith::topology::Jacobian::~Jacobian(void)
+{ // destructor
+ MatDestroy(*_matrix);
+ delete _matrix; _matrix = 0;
+} // destructor
+
+// ----------------------------------------------------------------------
+// Get PETSc matrix.
+const PetscMat*
+pylith::topology::Jacobian::matrix(void) const
+{ // matrix
+ return _matrix;
+} // matrix
+
+// ----------------------------------------------------------------------
+// Get PETSc matrix.
+PetscMat*
+pylith::topology::Jacobian::matrix(void)
+{ // matrix
+ return _matrix;
+} // matrix
+
+// ----------------------------------------------------------------------
+// Assemble matrix.
+void
+pylith::topology::Jacobian::assemble(const char* mode)
+{ // assemble
+ PetscErrorCode err = 0;
+ if (0 == strcmp(mode, "final_assembly")) {
+ err = MatAssemblyBegin(*_matrix, MAT_FINAL_ASSEMBLY);
+ if (err) {
+ PetscError(__LINE__,__FUNCT__,__FILE__,__SDIR__,err,0," ");
+ throw std::runtime_error("Error beginning final assembly of sparse "
+ "matrix associated with system Jacobian.");
+ } // if
+ err = MatAssemblyEnd(*_matrix, MAT_FINAL_ASSEMBLY);
+ if (err) {
+ PetscError(__LINE__,__FUNCT__,__FILE__,__SDIR__,err,0," ");
+ throw std::runtime_error("Error ending final assembly of sparse "
+ "matrix associated with system Jacobian.");
+ } // if
+ } else if (0 == strcmp(mode, "flush_assembly")) {
+ err = MatAssemblyBegin(*_matrix, MAT_FLUSH_ASSEMBLY);
+ if (err) {
+ PetscError(__LINE__,__FUNCT__,__FILE__,__SDIR__,err,0," ");
+ throw std::runtime_error("Error beginning flush assembly of sparse "
+ "matrix associated with system Jacobian.");
+ } // if
+ err = MatAssemblyEnd(*_matrix, MAT_FLUSH_ASSEMBLY);
+ if (err) {
+ PetscError(__LINE__,__FUNCT__,__FILE__,__SDIR__,err,0," ");
+ throw std::runtime_error("Error ending flush assembly of sparse "
+ "matrix associated with system Jacobian.");
+ } // if
+ } else
+ throw std::runtime_error("Unknown mode for assembly of sparse matrix "
+ "associated with system Jacobian.");
+} // assemble
+
+// ----------------------------------------------------------------------
+// Set entries in matrix to zero (retain structure).
+void
+pylith::topology::Jacobian::zero(void)
+{ // zero
+ PetscErrorCode err = MatZeroEntries(*_matrix);
+ if (err) {
+ PetscError(__LINE__,__FUNCT__,__FILE__,__SDIR__,err,0," ");
+ throw std::runtime_error("Error zeroing entries of sparse matrix "
+ "associated with system Jacobian.");
+ } // if
+} // zero
+
+// ----------------------------------------------------------------------
+// View matrix to stdout.
+void
+pylith::topology::Jacobian::view(void)
+{ // view
+ PetscErrorCode err = MatView(*_matrix, PETSC_VIEWER_STDOUT_WORLD);
+ if (err) {
+ PetscError(__LINE__,__FUNCT__,__FILE__,__SDIR__,err,0," ");
+ throw std::runtime_error("Error viewing sparse matrix associatd "
+ "with system Jacobian.");
+ } // if
+} // view
+
+// ----------------------------------------------------------------------
+// Write matrix to binary file.
+void
+pylith::topology::Jacobian::write(const char* filename)
+{ // write
+ PetscViewer viewer;
+
+ const MPI_Comm comm = _fields.mesh().comm();
+
+ PetscErrorCode err =
+ PetscViewerBinaryOpen(comm, filename,
+ FILE_MODE_WRITE, &viewer);
+ if (err) {
+ PetscError(__LINE__,__FUNCT__,__FILE__,__SDIR__,err,0," ");
+ throw std::runtime_error("Could not create PETSc binary viewer for "
+ "sparse matrix associated with system Jacobian.");
+ } // if
+
+ err = MatView(*_matrix, viewer);
+ if (err) {
+ PetscError(__LINE__,__FUNCT__,__FILE__,__SDIR__,err,0," ");
+ throw std::runtime_error("Could not view PETSc sparse matrix associated "
+ "with system Jacobian.");
+ } // if
+
+ err = PetscViewerDestroy(viewer);
+ if (err) {
+ PetscError(__LINE__,__FUNCT__,__FILE__,__SDIR__,err,0," ");
+ throw std::runtime_error("Could not destroy PETSc binary viewer for "
+ "sparse matrix associated with system Jacobian.");
+ } // if
+
+} // write
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/topology/Jacobian.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/topology/Jacobian.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/topology/Jacobian.hh 2009-03-12 23:10:01 UTC (rev 14310)
@@ -0,0 +1,91 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file libsrc/topology/Jacobian.hh
+ *
+ * @brief Jacobian of the system as a PETSc sparse matrix.
+ */
+
+#if !defined(pylith_topology_jacobian_hh)
+#define pylith_topology_jacobian_hh
+
+// Include directives ---------------------------------------------------
+#include "topologyfwd.hh" // forward declarations
+
+#include "pylith/utils/petscfwd.h" // HOLDSA PetscMat
+
+// Jacobian -------------------------------------------------------------
+class pylith::topology::Jacobian
+{ // Jacobian
+ friend class TestJacobian; // unit testing
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /** Default constructor.
+ *
+ * @param fields Fields associated with mesh and solution of the problem.
+ */
+ Jacobian(const SolutionFields& fields);
+
+ /// Destructor.
+ ~Jacobian(void);
+
+ /** Get PETSc matrix.
+ *
+ * @returns PETSc sparse matrix.
+ */
+ const PetscMat* matrix(void) const;
+
+ /** Get PETSc matrix.
+ *
+ * @returns PETSc sparse matrix.
+ */
+ PetscMat* matrix(void);
+
+ /** Assemble matrix.
+ *
+ * @param mode Assembly mode.
+ */
+ void assemble(const char* mode);
+
+ /// Set entries in matrix to zero (retain structure).
+ void zero(void);
+
+ /// View matrix to stdout.
+ void view(void);
+
+ /** Write matrix to binary file.
+ *
+ * @param filename Name of file.
+ */
+ void write(const char* filename);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ const SolutionFields& _fields; ///< Solution fields associated with problem.
+ PetscMat* _matrix; ///< Sparse matrix for Jacobian of problem.
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ Jacobian(const Jacobian&); ///< Not implemented
+ const Jacobian& operator=(const Jacobian&); ///< Not implemented
+
+}; // Jacobian
+
+#endif // pylith_topology_jacobian_hh
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/modulesrc/topology/Jacobian.i
===================================================================
--- short/3D/PyLith/branches/pylith-swig/modulesrc/topology/Jacobian.i (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/modulesrc/topology/Jacobian.i 2009-03-12 23:10:01 UTC (rev 14310)
@@ -0,0 +1,72 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file modulesrc/topology/Jacobian.i
+ *
+ * @brief Python interface to C++ Jacobian object.
+ */
+
+namespace pylith {
+ namespace topology {
+
+ class Jacobian
+ { // Jacobian
+
+ // PUBLIC MEMBERS /////////////////////////////////////////////////
+ public :
+
+ /** Default constructor.
+ *
+ * @param fields Fields associated with mesh and solution of the problem.
+ */
+ Jacobian(const SolutionFields& fields);
+
+ /// Destructor.
+ ~Jacobian(void);
+
+ /** Get PETSc matrix.
+ *
+ * @returns PETSc sparse matrix.
+ */
+ const PetscMat* matrix(void) const;
+
+ /** Get PETSc matrix.
+ *
+ * @returns PETSc sparse matrix.
+ */
+ PetscMat* matrix(void);
+
+ /** Assemble matrix.
+ *
+ * @param mode Assembly mode.
+ */
+ void assemble(const char* mode);
+
+ /// Set entries in matrix to zero (retain structure).
+ void zero(void);
+
+ /// View matrix to stdout.
+ void view(void);
+
+ /** Write matrix to binary file.
+ *
+ * @param filename Name of file.
+ */
+ void write(const char* filename);
+
+ }; // Jacobian
+
+ } // topology
+} // pylith
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/pylith/topology/Jacobian.py
===================================================================
--- short/3D/PyLith/branches/pylith-swig/pylith/topology/Jacobian.py (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/pylith/topology/Jacobian.py 2009-03-12 23:10:01 UTC (rev 14310)
@@ -0,0 +1,38 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/topology/Jacobian.py
+##
+## @brief Python object for system Jacobian.
+
+from topology import Jacobian as ModuleJacobian
+
+# ----------------------------------------------------------------------
+# Jacobian class
+class Jacobian(ModuleJacobian):
+ """
+ Python object for system Jacobian.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, fields):
+ """
+ Constructor.
+
+ @param fields Solution fields.
+ """
+ ModuleJacobian.__init__(self, fields)
+ return
+
+
+# End of file
Added: short/3D/PyLith/branches/pylith-swig/unittests/libtests/topology/TestJacobian.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/topology/TestJacobian.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/topology/TestJacobian.cc 2009-03-12 23:10:01 UTC (rev 14310)
@@ -0,0 +1,130 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+#include <portinfo>
+
+#include "TestJacobian.hh" // Implementation of class methods
+
+#include "pylith/topology/Jacobian.hh" // USES Jacobian
+
+#include "pylith/topology/Mesh.hh" // USES Mesh
+#include "pylith/topology/Field.hh" // USES Field
+#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
+
+#include "pylith/meshio/MeshIOAscii.hh" // USES MeshIOAscii
+
+// ----------------------------------------------------------------------
+CPPUNIT_TEST_SUITE_REGISTRATION( pylith::topology::TestJacobian );
+
+// ----------------------------------------------------------------------
+// Test constructor.
+void
+pylith::topology::TestJacobian::testConstructor(void)
+{ // testConstructor
+ Mesh mesh;
+ SolutionFields fields(mesh);
+ _initialize(&mesh, &fields);
+ Jacobian jacobian(fields);
+} // testConstructor
+
+// ----------------------------------------------------------------------
+// Test matrix().
+void
+pylith::topology::TestJacobian::testMatrix(void)
+{ // testMatrix
+ Mesh mesh;
+ SolutionFields fields(mesh);
+ _initialize(&mesh, &fields);
+ Jacobian jacobian(fields);
+
+ const PetscMat* matrix = jacobian.matrix();
+ CPPUNIT_ASSERT(0 != matrix);
+} // testMatrix
+
+// ----------------------------------------------------------------------
+// Test assemble().
+void
+pylith::topology::TestJacobian::testAssemble(void)
+{ // testAssemble
+ Mesh mesh;
+ SolutionFields fields(mesh);
+ _initialize(&mesh, &fields);
+ Jacobian jacobian(fields);
+
+ jacobian.assemble("flush_assembly");
+ jacobian.assemble("final_assembly");
+} // testAssemble
+
+// ----------------------------------------------------------------------
+// Test zero().
+void
+pylith::topology::TestJacobian::testZero(void)
+{ // testZero
+ Mesh mesh;
+ SolutionFields fields(mesh);
+ _initialize(&mesh, &fields);
+ Jacobian jacobian(fields);
+
+ jacobian.zero();
+} // testZero
+
+// ----------------------------------------------------------------------
+// Test view().
+void
+pylith::topology::TestJacobian::testView(void)
+{ // testView
+ Mesh mesh;
+ SolutionFields fields(mesh);
+ _initialize(&mesh, &fields);
+ Jacobian jacobian(fields);
+
+ jacobian.assemble("final_assembly");
+
+ jacobian.view();
+} // testView
+
+// ----------------------------------------------------------------------
+// Test write().
+void
+pylith::topology::TestJacobian::testWrite(void)
+{ // testWrite
+ Mesh mesh;
+ SolutionFields fields(mesh);
+ _initialize(&mesh, &fields);
+ Jacobian jacobian(fields);
+
+ jacobian.assemble("final_assembly");
+
+ jacobian.write("jacobian.mat");
+} // testWrite
+
+// ----------------------------------------------------------------------
+void
+pylith::topology::TestJacobian::_initialize(Mesh* mesh,
+ SolutionFields* fields) const
+{ // _initialize
+ CPPUNIT_ASSERT(0 != mesh);
+
+ meshio::MeshIOAscii iohandler;
+ iohandler.filename("data/tri3.mesh");
+ iohandler.read(mesh);
+
+ fields->add("disp t+dt");
+ fields->solutionName("disp t+dt");
+ Field<Mesh>& solution = fields->solution();
+ solution.newSection(FieldBase::VERTICES_FIELD, mesh->dimension());
+ solution.allocate();
+ solution.zero();
+} // _initialize
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/unittests/libtests/topology/TestJacobian.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/libtests/topology/TestJacobian.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/unittests/libtests/topology/TestJacobian.hh 2009-03-12 23:10:01 UTC (rev 14310)
@@ -0,0 +1,88 @@
+// -*- C++ -*-
+//
+// ----------------------------------------------------------------------
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ----------------------------------------------------------------------
+//
+
+/**
+ * @file unittests/libtests/topology/TestJacobian.hh
+ *
+ * @brief C++ TestJacobian object.
+ *
+ * C++ unit testing for Jacobian.
+ */
+
+#if !defined(pylith_topology_testjacobian_hh)
+#define pylith_topology_testjacobian_hh
+
+#include <cppunit/extensions/HelperMacros.h>
+
+#include "pylith/topology/topologyfwd.hh"
+
+/// Namespace for pylith package
+namespace pylith {
+ namespace topology {
+ class TestJacobian;
+ } // topology
+} // pylith
+
+/// C++ unit testing for Jacobian.
+class pylith::topology::TestJacobian : public CppUnit::TestFixture
+{ // class TestJacobian
+
+ // CPPUNIT TEST SUITE /////////////////////////////////////////////////
+ CPPUNIT_TEST_SUITE( TestJacobian );
+
+ CPPUNIT_TEST( testConstructor );
+ CPPUNIT_TEST( testMatrix );
+ CPPUNIT_TEST( testAssemble );
+ CPPUNIT_TEST( testZero );
+ CPPUNIT_TEST( testView );
+ CPPUNIT_TEST( testWrite );
+
+ CPPUNIT_TEST_SUITE_END();
+
+ // PUBLIC METHODS /////////////////////////////////////////////////////
+public :
+
+ /// Test constructor.
+ void testConstructor(void);
+
+ /// Test matrix().
+ void testMatrix(void);
+
+ /// Test assemble().
+ void testAssemble(void);
+
+ /// Test zero().
+ void testZero(void);
+
+ /// Test view().
+ void testView(void);
+
+ /// Test write().
+ void testWrite(void);
+
+ // PRIVATE METHODS ////////////////////////////////////////////////////
+private :
+
+ /** Initialize mesh for Jacobian.
+ *
+ * @param mesh Finite-element mesh.
+ * @param fields Solution fields.
+ */
+ void _initialize(Mesh* mesh,
+ SolutionFields* fields) const;
+
+}; // class TestJacobian
+
+#endif // pylith_topology_jacobian_hh
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/unittests/pytests/topology/TestJacobian.py
===================================================================
--- short/3D/PyLith/branches/pylith-swig/unittests/pytests/topology/TestJacobian.py (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/unittests/pytests/topology/TestJacobian.py 2009-03-12 23:10:01 UTC (rev 14310)
@@ -0,0 +1,137 @@
+#!/usr/bin/env python
+#
+# ======================================================================
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# {LicenseText}
+#
+# ======================================================================
+#
+
+## @file unittests/pytests/topology/TestJacobian.py
+
+## @brief Unit testing of Jacobian object.
+
+import unittest
+
+from pylith.topology.Jacobian import Jacobian
+
+
+# ----------------------------------------------------------------------
+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(normalizer, debug=False, interpolate=False)
+
+ from pylith.topology.SolutionFields import SolutionFields
+ fields = SolutionFields(self.mesh)
+ fields.add("disp t+dt")
+ 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(self.fields)
+ return
+
+
+ def test_constructor(self):
+ """
+ Test constructor.
+ """
+ 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.
+ """
+ self.jacobian.assemble("final_assembly")
+ self.jacobian.view()
+
+ # No testing 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.assemble("final_assembly")
+ self.jacobian.write("jacobian.mat")
+
+ # No testing of result.
+ return
+
+
+# End of file
More information about the CIG-COMMITS
mailing list