[cig-commits] r14261 - in short/3D/PyLith/branches/pylith-swig: . libsrc libsrc/problems
brad at geodynamics.org
brad at geodynamics.org
Sun Mar 8 13:02:11 PDT 2009
Author: brad
Date: 2009-03-08 13:02:11 -0700 (Sun, 08 Mar 2009)
New Revision: 14261
Added:
short/3D/PyLith/branches/pylith-swig/libsrc/problems/
short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.cc
short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.hh
short/3D/PyLith/branches/pylith-swig/libsrc/problems/Makefile.am
short/3D/PyLith/branches/pylith-swig/libsrc/problems/problemsfwd.hh
Modified:
short/3D/PyLith/branches/pylith-swig/configure.ac
short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
Log:
Added C++ Formulation object for interfacing with SNES.
Modified: short/3D/PyLith/branches/pylith-swig/configure.ac
===================================================================
--- short/3D/PyLith/branches/pylith-swig/configure.ac 2009-03-08 15:43:24 UTC (rev 14260)
+++ short/3D/PyLith/branches/pylith-swig/configure.ac 2009-03-08 20:02:11 UTC (rev 14261)
@@ -214,6 +214,7 @@
libsrc/faults/Makefile
libsrc/materials/Makefile
libsrc/meshio/Makefile
+ libsrc/problems/Makefile
libsrc/topology/Makefile
libsrc/utils/Makefile
modulesrc/Makefile
Modified: short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-03-08 15:43:24 UTC (rev 14260)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/Makefile.am 2009-03-08 20:02:11 UTC (rev 14261)
@@ -14,10 +14,11 @@
utils \
topology \
meshio \
+ materials
+ bc \
feassemble \
- bc \
faults \
- materials
+ problems
lib_LTLIBRARIES = libpylith.la
@@ -74,6 +75,7 @@
meshio/PsetFile.cc \
meshio/PsetFileAscii.cc \
meshio/PsetFileBinary.cc \
+ problems/Formulation.cc \
topology/FieldBase.cc \
topology/Mesh.cc \
topology/MeshOps.cc \
Added: short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.cc (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.cc 2009-03-08 20:02:11 UTC (rev 14261)
@@ -0,0 +1,177 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include <portinfo>
+
+#include "Formulation.hh" // implementation of class methods
+
+#include "pylith/topology/Mesh.hh" // USES Quadrature<Mesh>
+#include "pylith/topology/SubMesh.hh" // USES Quadrature<SubMesh>
+#include "pylith/feassemble/Quadrature.hh" // USES Integrator<Quadrature>
+#include "pylith/feassemble/Integrator.hh" // USES Integrator
+
+#include <cassert> // USES assert()
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::problems::Formulation::Formulation(void)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::problems::Formulation::~Formulation(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Generic C interface for reformResidual for integration with
+// PETSc SNES solvers.
+void
+pylith::problems::Formulation::reformResidual(void* context)
+{ // reformResidual
+ assert(0 != context);
+ ArgsResidual* args = (ArgsResidual*) context;
+ assert(0 != args);
+ assert(0 != args->object);
+ args->object->reformResidual(args->residual, args->fields, args->t, args->dt);
+} // reformResidual
+
+// ----------------------------------------------------------------------
+// Generic C interface for reformJacobian for integration with
+// PETSc SNES solvers.
+void
+pylith::problems::Formulation::reformJacobian(void* context)
+{ // reformJacobian
+ assert(0 != context);
+ ArgsJacobian* args = (ArgsJacobian*) context;
+ assert(0 != args);
+ assert(0 != args->object);
+ args->object->reformJacobian(args->jacobian, args->fields, args->t, args->dt);
+} // reformJacobian
+
+// ----------------------------------------------------------------------
+// Set integrators over the mesh.
+void
+pylith::problems::Formulation::meshIntegrators(IntegratorMesh** integrators,
+ const int numIntegrators)
+{ // meshIntegrators
+ assert( (0 == integrators && 0 == numIntegrators) ||
+ (0 != integrators && 0 < numIntegrators) );
+ _meshIntegrators.resize(numIntegrators);
+ for (int i=0; i < numIntegrators; ++i)
+ _meshIntegrators[i] = integrators[i];
+} // meshIntegrators
+
+// ----------------------------------------------------------------------
+// Set integrators over lower-dimension meshes.
+void
+pylith::problems::Formulation::submeshIntegrators(IntegratorSubMesh** integrators,
+ const int numIntegrators)
+{ // submeshIntegrators
+ assert( (0 == integrators && 0 == numIntegrators) ||
+ (0 != integrators && 0 < numIntegrators) );
+ _submeshIntegrators.resize(numIntegrators);
+ for (int i=0; i < numIntegrators; ++i)
+ _submeshIntegrators[i] = integrators[i];
+} // submeshIntegrators
+
+// ----------------------------------------------------------------------
+// Reform system residual.
+void
+pylith::problems::Formulation::reformResidual(
+ topology::Field<topology::Mesh>* const residual,
+ topology::SolutionFields* const fields,
+ const double t,
+ const double dt)
+{ // reformResidual
+ assert(0 != residual);
+ assert(0 != fields);
+
+ // Set residual to zero.
+ residual->zero();
+
+ // Add in contributions that require assembly.
+ int numIntegrators = _meshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i) {
+ _meshIntegrators[i]->timeStep(dt);
+ _meshIntegrators[i]->integrateResidual(*residual, t, fields);
+ } // for
+ numIntegrators = _submeshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i) {
+ _submeshIntegrators[i]->timeStep(dt);
+ _submeshIntegrators[i]->integrateResidual(*residual, t, fields);
+ } // for
+
+ // Assemble residual.
+ residual->complete();
+
+ // Add in contributions that do not require assembly.
+ numIntegrators = _meshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _meshIntegrators[i]->integrateResidualAssembled(*residual, t, fields);
+ numIntegrators = _submeshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _submeshIntegrators[i]->integrateResidual(*residual, t, fields);
+} // reformResidual
+
+// ----------------------------------------------------------------------
+// Reform system Jacobian.
+void
+pylith::problems::Formulation::reformJacobian(
+ PetscMat* jacobian,
+ topology::SolutionFields* const fields,
+ const double t,
+ const double dt)
+{ // reformJacobian
+ assert(0 != jacobian);
+ assert(0 != fields);
+
+ // Set residual to zero.
+ MatZeroEntries(*jacobian);
+
+ // Add in contributions that do not require assembly.
+ int numIntegrators = _meshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _meshIntegrators[i]->integrateJacobianAssembled(jacobian, t, fields);
+ numIntegrators = _submeshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _submeshIntegrators[i]->integrateJacobianAssembled(jacobian, t, fields);
+
+ // Assemble residual.
+ PetscErrorCode err = 0;
+ err = MatAssemblyBegin(*jacobian, MAT_FLUSH_ASSEMBLY);
+ if (0 != err)
+ throw std::runtime_error("Assembly of matrix failed.");
+ err = MatAssemblyEnd(*jacobian, MAT_FLUSH_ASSEMBLY);
+ if (0 != err)
+ throw std::runtime_error("Assembly of matrix failed.");
+
+ // Add in contributions that require assembly.
+ numIntegrators = _meshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _meshIntegrators[i]->integrateJacobian(jacobian, t, fields);
+ numIntegrators = _submeshIntegrators.size();
+ for (int i=0; i < numIntegrators; ++i)
+ _submeshIntegrators[i]->integrateJacobian(jacobian, t, fields);
+
+ // Assemble residual.
+ err = MatAssemblyBegin(*jacobian, MAT_FINAL_ASSEMBLY);
+ if (0 != err)
+ throw std::runtime_error("Assembly of matrix failed.");
+ err = MatAssemblyEnd(*jacobian, MAT_FINAL_ASSEMBLY);
+ if (0 != err)
+ throw std::runtime_error("Assembly of matrix failed.");
+} // reformJacobian
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.hh 2009-03-08 20:02:11 UTC (rev 14261)
@@ -0,0 +1,142 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/**
+ * @file pylith/problems/Formulation.hh
+ *
+ * @brief Object for managing forming Jacobian and residual for the problem.
+ */
+
+#if !defined(pylith_problems_formulation_hh)
+#define pylith_problems_formulation_hh
+
+// Include directives ---------------------------------------------------
+#include "problemsfwd.hh" // forward declarations
+
+#include "pylith/feassemble/feassemblefwd.hh" // USES Integrator
+#include "pylith/topology/topologyfwd.hh" // USES Mesh, Field, SolutionFields
+#include "pylith/utils/array.hh" // USES std::vector
+#include "pylith/utils/petscfwd.h" // USES PetscMat
+
+// Formulation ----------------------------------------------------------
+class pylith::problems::Formulation
+{ // Integrator
+ friend class TestFormulation; // unit testing
+
+// PUBLIC STRUCTS ///////////////////////////////////////////////////////
+public :
+
+ struct ArgsResidual {
+ Formulation* object;
+ topology::Field<topology::Mesh>* const residual;
+ topology::SolutionFields* const fields;
+ double t;
+ double dt;
+ }; // ArgsResidual
+
+ struct ArgsJacobian {
+ Formulation* object;
+ PetscMat* jacobian;
+ topology::SolutionFields* const fields;
+ double t;
+ double dt;
+ }; // ArgsJacobian
+
+// PRIVATE TYPEDEFS /////////////////////////////////////////////////////
+private :
+
+ typedef feassemble::Integrator<feassemble::Quadrature<topology::Mesh> > IntegratorMesh;
+ typedef feassemble::Integrator<feassemble::Quadrature<topology::SubMesh> > IntegratorSubMesh;
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ Formulation(void);
+
+ /// Destructor
+ ~Formulation(void);
+
+ /** Generic C interface for reformResidual for integration with
+ * PETSc SNES solvers.
+ *
+ * @param context ArgsResidual structure with arguments.
+ */
+ static
+ void reformResidual(void* context);
+
+ /** Generic C interface for reformJacobian for integration with
+ * PETSc SNES solvers.
+ *
+ * @param context ArgsJacobian structure with arguments.
+ */
+ static
+ void reformJacobian(void* context);
+
+ /** Set integrators over the mesh.
+ *
+ * @param integrators Integrators over the mesh.
+ * @param numIntegrators Number of integrators.
+ */
+ void meshIntegrators(IntegratorMesh** integrators,
+ const int numIntegrators);
+
+ /** Set integrators over lower-dimension meshes.
+ *
+ * @param integrators Integrators over lower-dimension meshes.
+ * @param numIntegrators Number of integrators.
+ */
+ void submeshIntegrators(IntegratorSubMesh** integrators,
+ const int numIntegrators);
+
+ /** Reform system residual.
+ *
+ * @param residual Field containing values for residual
+ * @param fields Solution fields
+ * @param t Current time.
+ * @param dt Current time step (t -> t+dt).
+ */
+ void reformResidual(topology::Field<topology::Mesh>* const residual,
+ topology::SolutionFields* const fields,
+ const double t,
+ const double dt);
+
+ /** Reform system Jacobian.
+ *
+ * @param jacobian Sparse matrix for Jacobian of system.
+ * @param fields Solution fields
+ * @param t Current time.
+ * @param dt Current time step (t -> t+dt).
+ */
+ void reformJacobian(PetscMat* jacobian,
+ topology::SolutionFields* const fields,
+ const double t,
+ const double dt);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ std::vector<IntegratorMesh*> _meshIntegrators;
+ std::vector<IntegratorSubMesh*> _submeshIntegrators;
+
+// NOT IMPLEMENTED //////////////////////////////////////////////////////
+private :
+
+ Formulation(const Formulation&); ///< Not implemented
+ const Formulation& operator=(const Formulation&); ///< Not implemented
+
+}; // Formulation
+
+#endif // pylith_problems_formulation_hh
+
+
+// End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/problems/Makefile.am
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/problems/Makefile.am (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/problems/Makefile.am 2009-03-08 20:02:11 UTC (rev 14261)
@@ -0,0 +1,29 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+subpackage = problems
+include $(top_srcdir)/subpackage.am
+
+subpkginclude_HEADERS = \
+ Formulation.hh \
+ Formulation.icc \
+ problemsfwd.hh
+
+noinst_HEADERS =
+
+# export
+clean-local: clean-subpkgincludeHEADERS
+BUILT_SOURCES = export-subpkgincludeHEADERS
+CLEANFILES = export-subpkgincludeHEADERS
+
+
+# End of file
Added: short/3D/PyLith/branches/pylith-swig/libsrc/problems/problemsfwd.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/problems/problemsfwd.hh (rev 0)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/problems/problemsfwd.hh 2009-03-08 20:02:11 UTC (rev 14261)
@@ -0,0 +1,36 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+/** @file libsrc/problems/problemsfwd.hh
+ *
+ * @brief Forward declarations for PyLith problems objects.
+ *
+ * Including this header file eliminates the need to use separate
+ * forward declarations.
+ */
+
+#if !defined(pylith_problems_problemsfwd_hh)
+#define pylith_problems_problemsfwd_hh
+
+namespace pylith {
+ namespace problems {
+
+ class Formulation;
+
+ } // problems
+} // pylith
+
+
+#endif // pylith_problems_problemsfwd_hh
+
+
+// End of file
More information about the CIG-COMMITS
mailing list