[cig-commits] r14347 - short/3D/PyLith/branches/pylith-swig/libsrc/problems

[brad at geodynamics.org]

Author: brad
Date: 2009-03-16 15:10:28 -0700 (Mon, 16 Mar 2009)
New Revision: 14347

Modified:
   short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.cc
   short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.hh
Log:
More work on setting up Formulation.

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.cc
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.cc	2009-03-16 22:10:06 UTC (rev 14346)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.cc	2009-03-16 22:10:28 UTC (rev 14347)
@@ -19,7 +19,9 @@
 #include "pylith/feassemble/Quadrature.hh" // USES Integrator<Quadrature>
 #include "pylith/feassemble/Integrator.hh" // USES Integrator
 #include "pylith/topology/Jacobian.hh" // USES Jacobian
+#include "pylith/topology/SolutionFields.hh" // USES SolutionFields
 
+#include "pylith/utils/petscerror.h" // USES CHECK_PETSC_ERROR
 #include <cassert> // USES assert()
 
 // ----------------------------------------------------------------------
@@ -38,20 +40,67 @@
 // Generic C interface for reformResidual for integration with
 // PETSc SNES solvers.
 void
-pylith::problems::Formulation::reformResidual(void* context)
+pylith::problems::Formulation::reformResidual(PetscSNES snes,
+					      PetscVec solutionVec,
+					      PetscVec residualVec,
+					      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);  
+  assert(0 != args->fields);
+
+  PetscVec localVec;
+  PetscErrorCode err = 0;
+
+#if 0
+  PetscVecScatter scatter = fields->scatter();
+#else
+  VecScatter scatter = 0;
+#endif
+
+  // Copy solution information from PETSc vector into field
+  const ALE::Obj<topology::Mesh::RealSection>& solutionSection = 
+    args->fields->solution().section();
+  err = VecCreateSeqWithArray(PETSC_COMM_SELF,
+			      solutionSection->sizeWithBC(),
+			      solutionSection->restrictSpace(),
+			      &localVec); CHECK_PETSC_ERROR(err);
+  err = VecScatterBegin(scatter, solutionVec, localVec, INSERT_VALUES,
+			SCATTER_REVERSE); CHECK_PETSC_ERROR(err);
+  err = VecScatterEnd(scatter, solutionVec, localVec, INSERT_VALUES,
+		      SCATTER_REVERSE); CHECK_PETSC_ERROR(err);
+  err = VecDestroy(localVec); CHECK_PETSC_ERROR(err);
+
+  // Reform residual
+  args->object->reformResidual(args->residual, args->fields, 
+			       args->t, args->dt);  
+
+  // Copy residual information from field into PETSc vector
+  const ALE::Obj<topology::Mesh::RealSection>& residualSection = 
+    args->fields->get("residual").section();
+  err = VecCreateSeqWithArray(PETSC_COMM_SELF, 
+			      residualSection->sizeWithBC(),
+			      residualSection->restrictSpace(), 
+			      &localVec); CHECK_PETSC_ERROR(err);
+  err = VecScatterBegin(scatter, localVec, residualVec, INSERT_VALUES,
+			SCATTER_FORWARD); CHECK_PETSC_ERROR(err);
+  err = VecScatterEnd(scatter, localVec, residualVec, INSERT_VALUES,
+		      SCATTER_FORWARD); CHECK_PETSC_ERROR(err);
+  err = VecDestroy(localVec); CHECK_PETSC_ERROR(err);
 } // reformResidual
 
 // ----------------------------------------------------------------------
 // Generic C interface for reformJacobian for integration with
 // PETSc SNES solvers.
 void
-pylith::problems::Formulation::reformJacobian(void* context)
+pylith::problems::Formulation::reformJacobian(PetscSNES snes,
+					      PetscVec solutionVec,
+					      PetscMat jacobianMat,
+					      PetscMat preconditionerMat,
+					      int* preconditionerLayout,
+					      void* context)
 { // reformJacobian
   assert(0 != context);
   ArgsJacobian* args = (ArgsJacobian*) context;

Modified: short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.hh
===================================================================
--- short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.hh	2009-03-16 22:10:06 UTC (rev 14346)
+++ short/3D/PyLith/branches/pylith-swig/libsrc/problems/Formulation.hh	2009-03-16 22:10:28 UTC (rev 14347)
@@ -24,8 +24,12 @@
 
 #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 PetscVec, PetscMat, PetscSNES
+
+#include "pylith/utils/array.hh" // HASA std::vector
+
+
 // Formulation ----------------------------------------------------------
 class pylith::problems::Formulation
 { // Integrator
@@ -68,18 +72,34 @@
   /** Generic C interface for reformResidual for integration with
    * PETSc SNES solvers.
    *
+   * @param snes PETSc scalable nonlinear equation solver.
+   * @param solutionVec PETSc vector for solution.
+   * @param residualVec PETSc vector for residual.
    * @param context ArgsResidual structure with arguments.
    */
   static
-  void reformResidual(void* context);
+  void reformResidual(PetscSNES snes,
+		      PetscVec solutionVec,
+		      PetscVec residualVec,
+		      void* context);
 
   /** Generic C interface for reformJacobian for integration with
    * PETSc SNES solvers.
    *
+   * @param snes PETSc scalable nonlinear equation solver.
+   * @param solutionVec PETSc vector for solution.
+   * @param jacobianMat PETSc sparse matrix for system Jacobian.
+   * @param preconditionerMat PETSc sparse matrix for preconditioner.
+   * @param Flag indicating layout of preconditioner matrix.
    * @param context ArgsJacobian structure with arguments.
    */
   static
-  void reformJacobian(void* context);
+  void reformJacobian(PetscSNES snes,
+		      PetscVec solutionVec,
+		      PetscMat jacobianMat,
+		      PetscMat preconditionerMat,
+		      int* preconditionerLayout,
+		      void* context);
 
   /** Set integrators over the mesh.
    *