Commits

Jed Brown committed 97ea565 Merge with conflicts

Merge branch 'tautges/dmmoab'

* 'tautges/dmmoab'
DMMoab: Remove OBJSC from makefile
TS ex30.cxx: write output file matching example number (ex30-final.h5m)
DMMoab: alltest run TS ex30
DMMoab: ops pointers initialized in DMCreate_Moab
DMMoab: fix test example makefile
DMMoab: close memory leaks by explicitly calling C++ destructors
DMMoab: interface does not require that all of PETSc be built using C++
DMMoab: include <sstream> to use ostringstream
DMMoab: do not mix CHKERRXX() with CHKERRQ()
Fixing uninitialized memory, and adding output.
Fixing compile error in dmmoab.
Changes suggested by Jed and Barry.
Changes suggested by Matt in pull request.
Adding DMMoab and VecMoab, along with various examples and tests.

Conflicts:
src/ts/examples/tutorials/makefile

  • Participants
  • Parent commits 668511f, 6ac2e81

Comments (0)

Files changed (15)

 testexamples_CUSP: ${TESTEXAMPLES_CUSP}
 testexamples_YAML: ${TESTEXAMPLES_YAML}
 testexamples_THREADCOMM: ${TESTEXAMPLES_THREADCOMM}
+testexamples_MOAB: ${TESTEXAMPLES_MOAB}
 testexamples_X:
 testexamples_OPENGL:
 testexamples_MPE:
 	-@${OMAKE} testexamples_YAML TESTEXAMPLES_YAML=`echo ${TESTEXAMPLES_YAML} | sed s/runex[0-9]*[a-z0-9_]*//g`
 buildexamples_THREADCOMM:
 	-@${OMAKE} testexamples_THREADCOMM TESTEXAMPLES_THREADCOMM=`echo ${TESTEXAMPLES_THREADCOMM} | sed s/runex[0-9]*[a-z0-9_]*//g`
+buildexamples_MOAB:
+	-@${OMAKE} testexamples_MOAB TESTEXAMPLES_MOAB=`echo ${TESTEXAMPLES_MOAB} | sed s/runex[0-9]*[a-z0-9_]*//g`
 buildexamples_X:
 buildexamples_OPENGL:
 buildexamples_MPE:

File include/petscdm.h

 #define DMREDUNDANT "redundant"
 #define DMAKKT      "akkt"
 #define DMPATCH     "patch"
+#define DMMOAB      "moab"
 
 PETSC_EXTERN PetscFunctionList DMList;
 PETSC_EXTERN PetscBool         DMRegisterAllCalled;

File include/petscdmmoab.h

+#if !defined(__PETSCMOAB_H)
+#define __PETSCMOAB_H
+
+#include <petscvec.h>
+#include <petscmat.h>
+#include <petscdm.h>
+
+#include <moab/Core.hpp>
+#include <moab/ParallelComm.hpp>
+
+/* The MBERR macro is used to save typing. It checks a MOAB error code
+ * (rval) and calls SETERRQ if not MB_SUCCESS. A message (msg) can
+ * also be passed in. */
+#define MBERR(msg,rval) do{if(rval != moab::MB_SUCCESS) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"MOAB ERROR (%i): %s",rval,msg);} while(0)
+#define MBERRNM(rval) do{if(rval != moab::MB_SUCCESS) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"MOAB ERROR (%i)",rval);} while(0)
+
+PETSC_EXTERN PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *moab);
+PETSC_EXTERN PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::ParallelComm *pcomm, moab::Tag ltog_tag, moab::Range *range, DM *moab);
+PETSC_EXTERN PetscErrorCode DMMoabSetParallelComm(DM dm,moab::ParallelComm *pcomm);
+PETSC_EXTERN PetscErrorCode DMMoabGetParallelComm(DM dm,moab::ParallelComm **pcomm);
+PETSC_EXTERN PetscErrorCode DMMoabSetInterface(DM dm,moab::Interface *iface);
+PETSC_EXTERN PetscErrorCode DMMoabGetInterface(DM dm,moab::Interface **mbint);
+PETSC_EXTERN PetscErrorCode DMMoabSetRange(DM dm,moab::Range range);
+PETSC_EXTERN PetscErrorCode DMMoabGetRange(DM dm,moab::Range *range);
+PETSC_EXTERN PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm,moab::Tag ltogtag);
+PETSC_EXTERN PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm,moab::Tag *ltog_tag);
+
+PETSC_EXTERN PetscErrorCode DMMoabCreateVector(DM dm, moab::Tag tag,PetscInt tag_size,moab::Range range,PetscBool serial, PetscBool destroy_tag,Vec *X);
+PETSC_EXTERN PetscErrorCode DMMoabGetVecTag(Vec vec,moab::Tag *tag);
+PETSC_EXTERN PetscErrorCode DMMoabGetVecRange(Vec vec,moab::Range *range);
+#endif

File src/dm/impls/makefile

 ALL: lib
 
-DIRS     = da adda sliced composite mesh cartesian redundant plex shell patch
+DIRS     = da adda sliced composite mesh cartesian redundant plex shell patch moab
 LOCDIR   = src/dm/impls
 MANSEC   = DM
 

File src/dm/impls/moab/dmmoab.cxx

+#include <petsc-private/dmimpl.h> /*I  "petscdm.h"   I*/
+#include <petsc-private/vecimpl.h> /*I  "petscdm.h"   I*/
+
+#include <petscdmmoab.h>
+#include <MBTagConventions.hpp>
+#include <sstream>
+
+typedef struct {
+  PetscInt bs; /* Number of degrees of freedom on each entity, aka tag size in moab */
+  PetscBool icreatedinstance; /* true if DM created moab instance internally, will destroy instance in DMDestroy */
+  moab::ParallelComm *pcomm;
+  moab::Interface *mbiface;
+  moab::Tag ltog_tag; /* moab supports "global id" tags, which are usually local to global numbering */
+  moab::Range range;
+} DM_Moab;
+
+typedef struct {
+  moab::Interface    *mbiface;
+  moab::ParallelComm *pcomm;
+  moab::Range         tag_range; /* entities to which this tag applies */
+  moab::Tag           tag;
+  moab::Tag           ltog_tag;
+  PetscInt            tag_size;
+  PetscBool           new_tag;
+  PetscBool           serial;
+
+} Vec_MOAB;
+
+#undef __FUNCT__
+#define __FUNCT__ "DMCreateGlobalVector_Moab"
+PetscErrorCode DMCreateGlobalVector_Moab(DM dm,Vec *gvec)
+{
+  PetscErrorCode  ierr;
+  DM_Moab         *dmmoab = (DM_Moab*)dm->data;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  PetscValidPointer(gvec,2);
+  PetscInt block_size = ((DM_Moab*)dm->data)->bs;
+  moab::Tag tag = 0;
+  ierr = DMMoabCreateVector(dm,tag,block_size,dmmoab->range,PETSC_FALSE,PETSC_TRUE,gvec);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMCreateLocalVector_Moab"
+PetscErrorCode DMCreateLocalVector_Moab(DM dm,Vec *gvec)
+{
+  PetscErrorCode  ierr;
+  DM_Moab         *dmmoab = (DM_Moab*)dm->data;
+
+  PetscFunctionBegin;
+  PetscInt bs = 1;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  PetscValidPointer(gvec,2);
+  moab::Tag tag = 0;
+  ierr = DMMoabCreateVector(dm,tag,bs,dmmoab->range,PETSC_TRUE,PETSC_TRUE,gvec);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMDestroy_Moab"
+PetscErrorCode DMDestroy_Moab(DM dm)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  if (((DM_Moab*)dm->data)->icreatedinstance) {
+    delete ((DM_Moab*)dm->data)->mbiface;
+    ((DM_Moab*)dm->data)->mbiface = NULL;
+    ((DM_Moab*)dm->data)->pcomm = NULL;
+    ((DM_Moab*)dm->data)->range.~Range();
+  }
+  ierr = PetscFree(dm->data);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMCreate_Moab"
+PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
+{
+  DM_Moab        *moab;
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  ierr = PetscNewLog(dm,DM_Moab,&moab);CHKERRQ(ierr);
+  dm->data = moab;
+  new (moab) DM_Moab();
+
+  dm->ops->createglobalvector              = DMCreateGlobalVector_Moab;
+  dm->ops->createlocalvector               = DMCreateLocalVector_Moab;
+  dm->ops->destroy                         = DMDestroy_Moab;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabCreate"
+/*@
+  DMMoabCreate - Creates a DMMoab object, which encapsulates a moab instance
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. comm - The communicator for the DMMoab object
+
+  Output Parameter:
+. moab  - The DMMoab object
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *moab)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  PetscValidPointer(moab,2);
+  ierr = DMCreate(comm, moab);CHKERRQ(ierr);
+  ierr = DMSetType(*moab, DMMOAB);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabCreateMoab"
+/*@
+  DMMoabCreate - Creates a DMMoab object, optionally from an instance and other data
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. comm - The communicator for the DMMoab object
+. moab - (ptr to) the MOAB Instance; if passed in NULL, MOAB instance is created inside PETSc, and destroyed
+         along with the DMMoab
+. pcomm - (ptr to) a ParallelComm; if NULL, creates one internally for the whole communicator
+. ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
+. range - If non-NULL, contains range of entities to which DOFs will be assigned
+
+  Output Parameter:
+. moab  - The DMMoab object
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::ParallelComm *pcomm, moab::Tag ltog_tag, moab::Range *range, DM *moab)
+{
+  PetscErrorCode ierr;
+  DM_Moab        *dmmoab;
+
+  PetscFunctionBegin;
+  PetscValidPointer(moab,2);
+  ierr = DMMoabCreate(comm, moab);CHKERRQ(ierr);
+  dmmoab = (DM_Moab*)(*moab)->data;
+
+  if (!mbiface) {
+    mbiface = new moab::Core();
+    dmmoab->icreatedinstance = PETSC_TRUE;
+  }
+  else
+    dmmoab->icreatedinstance = PETSC_FALSE;
+
+  if (!pcomm) {
+    PetscInt rank, nprocs;
+    MPI_Comm_rank(comm, &rank);
+    MPI_Comm_size(comm, &nprocs);
+    pcomm = new moab::ParallelComm(mbiface, comm);
+  }
+
+    // do the initialization of the DM
+  dmmoab->bs       = 0;
+  dmmoab->pcomm    = pcomm;
+  dmmoab->mbiface    = mbiface;
+  dmmoab->ltog_tag = ltog_tag;
+
+  ierr = DMMoabSetInterface(*moab, mbiface);CHKERRQ(ierr);
+  if (!pcomm) pcomm = new moab::ParallelComm(mbiface, comm);
+  ierr = DMMoabSetParallelComm(*moab, pcomm);CHKERRQ(ierr);
+  if (!ltog_tag) {
+    moab::ErrorCode merr = mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, ltog_tag);MBERRNM(merr);
+  }
+  if (ltog_tag) {
+    ierr = DMMoabSetLocalToGlobalTag(*moab, ltog_tag);CHKERRQ(ierr);
+  }
+  if (range) {
+    ierr = DMMoabSetRange(*moab, *range);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabSetParallelComm"
+/*@
+  DMMoabSetParallelComm - Set the ParallelComm used with this DMMoab
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. dm    - The DMMoab object being set
+. pcomm - The ParallelComm being set on the DMMoab
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabSetParallelComm(DM dm,moab::ParallelComm *pcomm)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  ((DM_Moab*)dm->data)->pcomm = pcomm;
+  ((DM_Moab*)dm->data)->mbiface = pcomm->get_moab();
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabGetParallelComm"
+/*@
+  DMMoabGetParallelComm - Get the ParallelComm used with this DMMoab
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. dm    - The DMMoab object being set
+
+  Output Parameter:
+. pcomm - The ParallelComm for the DMMoab
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabGetParallelComm(DM dm,moab::ParallelComm **pcomm)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  *pcomm = ((DM_Moab*)dm->data)->pcomm;
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabSetInterface"
+/*@
+  DMMoabSetInterface - Set the MOAB instance used with this DMMoab
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. dm      - The DMMoab object being set
+. mbiface - The MOAB instance being set on this DMMoab
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabSetInterface(DM dm,moab::Interface *mbiface)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  ((DM_Moab*)dm->data)->pcomm = NULL;
+  ((DM_Moab*)dm->data)->mbiface = mbiface;
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabGetInterface"
+/*@
+  DMMoabGetInterface - Get the MOAB instance used with this DMMoab
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. dm      - The DMMoab object being set
+
+  Output Parameter:
+. mbiface - The MOAB instance set on this DMMoab
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabGetInterface(DM dm,moab::Interface **mbiface)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  *mbiface = ((DM_Moab*)dm->data)->mbiface;
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabSetRange"
+/*@
+  DMMoabSetRange - Set the entities having DOFs on this DMMoab
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. dm    - The DMMoab object being set
+. range - The entities treated by this DMMoab
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabSetRange(DM dm,moab::Range range)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  ((DM_Moab*)dm->data)->range = range;
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabGetRange"
+/*@
+  DMMoabGetRange - Get the entities having DOFs on this DMMoab
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. dm    - The DMMoab object being set
+
+  Output Parameter:
+. range - The entities treated by this DMMoab
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabGetRange(DM dm,moab::Range *range)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  *range = ((DM_Moab*)dm->data)->range;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabSetLocalToGlobalTag"
+/*@
+  DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. dm      - The DMMoab object being set
+. ltogtag - The MOAB tag used for local to global ids
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm,moab::Tag ltogtag)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  ((DM_Moab*)dm->data)->ltog_tag = ltogtag;
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabGetLocalToGlobalTag"
+/*@
+  DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. dm      - The DMMoab object being set
+
+  Output Parameter:
+. ltogtag - The MOAB tag used for local to global ids
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm,moab::Tag *ltog_tag)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  *ltog_tag = ((DM_Moab*)dm->data)->ltog_tag;
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabSetBlockSize"
+/*@
+  DMMoabSetBlockSize - Set the block size used with this DMMoab
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. dm - The DMMoab object being set
+. bs - The block size used with this DMMoab
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabSetBlockSize(DM dm,PetscInt bs)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  ((DM_Moab*)dm->data)->bs = bs;
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabGetBlockSize"
+/*@
+  DMMoabGetBlockSize - Get the block size used with this DMMoab
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. dm - The DMMoab object being set
+
+  Output Parameter:
+. bs - The block size used with this DMMoab
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabGetBlockSize(DM dm,PetscInt *bs)
+{
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(dm,DM_CLASSID,1);
+  *bs = ((DM_Moab*)dm->data)->bs;
+  PetscFunctionReturn(0);
+}
+
+
+// declare for use later but before they're defined
+PetscErrorCode DMMoab_VecUserDestroy(void *user);
+PetscErrorCode DMMoab_VecDuplicate(Vec x,Vec *y);
+PetscErrorCode DMMoab_CreateTagName(const moab::ParallelComm *pcomm,std::string& tag_name);
+PetscErrorCode DMMoab_CreateVector(moab::Interface *iface,moab::ParallelComm *pcomm,moab::Tag tag,PetscInt tag_size,moab::Tag ltog_tag,moab::Range range,PetscBool serial, PetscBool destroy_tag,Vec *vec);
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabCreateVector"
+/*@
+  DMMoabCreateVector - Create a Vec from either an existing tag, or a specified tag size, and a range of entities
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. dm          - The DMMoab object being set
+. tag         - If non-zero, block size will be taken from the tag size
+. tag_size    - If tag was zero, this parameter specifies the block size; unique tag name will be generated automatically
+. range       - If non-empty, Vec corresponds to these entities, otherwise to the entities set on the DMMoab
+. serial      - If true, this is a serial Vec, otherwise a parallel one
+. destroy_tag - If true, MOAB tag is destroyed with Vec, otherwise it is left on MOAB
+
+  Output Parameter:
+. vec         - The created vector
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabCreateVector(DM dm,moab::Tag tag,PetscInt tag_size,moab::Range range,PetscBool serial, PetscBool destroy_tag,Vec *vec)
+{
+  PetscErrorCode     ierr;
+
+  PetscFunctionBegin;
+
+  DM_Moab *dmmoab = (DM_Moab*)dm->data;
+  moab::ParallelComm *pcomm = dmmoab->pcomm;
+  moab::Interface *mbiface = dmmoab->mbiface;
+  moab::Tag ltog_tag = dmmoab->ltog_tag;
+
+  if (!tag && !tag_size) {
+    PetscFunctionReturn(PETSC_ERR_ARG_WRONG);
+  }
+  else {
+    ierr = DMMoab_CreateVector(mbiface,pcomm,tag,tag_size,ltog_tag,range,serial,destroy_tag,vec);CHKERRQ(ierr);
+  }
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoab_CreateVector"
+PetscErrorCode DMMoab_CreateVector(moab::Interface *mbiface,moab::ParallelComm *pcomm,moab::Tag tag,PetscInt tag_size,moab::Tag ltog_tag,moab::Range range,PetscBool serial, PetscBool destroy_tag,Vec *vec)
+{
+  PetscErrorCode     ierr;
+  moab::ErrorCode    merr;
+
+  PetscFunctionBegin;
+
+  if (!tag) {
+    std::string tag_name;
+    ierr = DMMoab_CreateTagName(pcomm,tag_name);CHKERRQ(ierr);
+
+      // Create the default value for the tag (all zeros):
+    std::vector<PetscScalar> default_value(tag_size, 0.0);
+
+      // Create the tag:
+    merr = mbiface->tag_get_handle(tag_name.c_str(),tag_size,moab::MB_TYPE_DOUBLE,tag,
+                                   moab::MB_TAG_DENSE | moab::MB_TAG_CREAT,default_value.data());MBERRNM(merr);
+  }
+  else {
+
+      // Make sure the tag data is of type "double":
+    moab::DataType tag_type;
+    merr = mbiface->tag_get_data_type(tag, tag_type);MBERRNM(merr);
+    if(tag_type != moab::MB_TYPE_DOUBLE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Tag data type must be MB_TYPE_DOUBLE");
+  }
+
+    // Create the MOAB internal data object
+  Vec_MOAB *vmoab;
+  ierr = PetscMalloc(sizeof(Vec_MOAB),&vmoab);CHKERRQ(ierr);
+  new (vmoab) Vec_MOAB();
+  vmoab->tag = tag;
+  vmoab->ltog_tag = ltog_tag;
+  vmoab->mbiface = mbiface;
+  vmoab->pcomm = pcomm;
+  vmoab->tag_range = range;
+  vmoab->new_tag = destroy_tag;
+  vmoab->serial = serial;
+  merr = mbiface->tag_get_length(tag,vmoab->tag_size);MBERR("tag_get_size", merr);
+
+    // Call tag_iterate. This will cause MOAB to allocate memory for the
+    // tag data if it hasn't already happened:
+  int  count;
+  void *void_ptr;
+  merr = mbiface->tag_iterate(tag,range.begin(),range.end(),count,void_ptr);MBERRNM(merr);
+
+    // Check to make sure the tag data is in a single sequence:
+  if ((unsigned)count != range.size()) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only create MOAB Vector for single sequence");
+  PetscScalar *data_ptr = (PetscScalar*)void_ptr;
+
+    // Create the PETSc Vector:
+  if(!serial) {
+      // This is an MPI Vector:
+    ierr = VecCreateMPIWithArray(vmoab->pcomm->comm(),vmoab->tag_size,vmoab->tag_size*range.size(),
+                                 PETSC_DECIDE,data_ptr,vec);CHKERRQ(ierr);
+
+      // Vector created, manually set local to global mapping:
+    ISLocalToGlobalMapping ltog;
+    PetscInt               *gindices = new PetscInt[range.size()];
+    PetscInt               count = 0;
+    moab::Range::iterator  iter;
+    for(iter = range.begin(); iter != range.end(); iter++) {
+      int dof;
+      merr = mbiface->tag_get_data(ltog_tag,&(*iter),1,&dof);MBERRNM(merr);
+      gindices[count] = dof;
+      count++;
+    }
+
+    ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,range.size(),gindices,
+                                        PETSC_COPY_VALUES,&ltog);CHKERRQ(ierr);
+    ierr = VecSetLocalToGlobalMappingBlock(*vec,ltog);CHKERRQ(ierr);
+
+      // Clean up:
+    ierr = ISLocalToGlobalMappingDestroy(&ltog);CHKERRQ(ierr);
+    delete [] gindices;
+  } else {
+      // This is a serial vector:
+    ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,vmoab->tag_size,vmoab->tag_size*range.size(),data_ptr,vec);CHKERRQ(ierr);
+  }
+
+
+  PetscContainer moabdata;
+  ierr = PetscContainerCreate(PETSC_COMM_SELF,&moabdata);CHKERRQ(ierr);
+  ierr = PetscContainerSetPointer(moabdata,vmoab);CHKERRQ(ierr);
+  ierr = PetscContainerSetUserDestroy(moabdata,DMMoab_VecUserDestroy);CHKERRQ(ierr);
+  ierr = PetscObjectCompose((PetscObject)*vec,"MOABData",(PetscObject)moabdata);CHKERRQ(ierr);
+  (*vec)->ops->duplicate = DMMoab_VecDuplicate;
+
+  ierr = PetscContainerDestroy(&moabdata);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabGetVecTag"
+/*@
+  DMMoabGetVecTag - Get the MOAB tag associated with this Vec
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. vec - Vec being queried
+
+  Output Parameter:
+. tag - Tag associated with this Vec
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabGetVecTag(Vec vec,moab::Tag *tag)
+{
+  PetscContainer  moabdata;
+  Vec_MOAB        *vmoab;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBegin;
+
+  // Get the MOAB private data:
+  ierr = PetscObjectQuery((PetscObject)vec,"MOABData", (PetscObject*) &moabdata);CHKERRQ(ierr);
+  ierr = PetscContainerGetPointer(moabdata, (void**) &vmoab);CHKERRQ(ierr);
+
+  *tag = vmoab->tag;
+
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoabGetVecRange"
+/*@
+  DMMoabGetVecRange - Get the MOAB entities associated with this Vec
+
+  Collective on MPI_Comm
+
+  Input Parameter:
+. vec   - Vec being queried
+
+  Output Parameter:
+. range - Entities associated with this Vec
+
+  Level: beginner
+
+.keywords: DMMoab, create
+@*/
+PetscErrorCode DMMoabGetVecRange(Vec vec,moab::Range *range)
+{
+  PetscContainer  moabdata;
+  Vec_MOAB        *vmoab;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBegin;
+
+  // Get the MOAB private data:
+  ierr = PetscObjectQuery((PetscObject)vec,"MOABData", (PetscObject*) &moabdata);CHKERRQ(ierr);
+  ierr = PetscContainerGetPointer(moabdata, (void**) &vmoab);CHKERRQ(ierr);
+
+  *range = vmoab->tag_range;
+
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoab_VecDuplicate"
+PetscErrorCode DMMoab_VecDuplicate(Vec x,Vec *y)
+{
+  PetscErrorCode ierr;
+  PetscFunctionBegin;
+  PetscValidHeaderSpecific(x,VEC_CLASSID,1);
+  PetscValidPointer(y,2);
+
+  // Get the Vec_MOAB struct for the original vector:
+  PetscContainer  moabdata;
+  Vec_MOAB        *vmoab;
+  ierr = PetscObjectQuery((PetscObject)x,"MOABData", (PetscObject*) &moabdata);CHKERRQ(ierr);
+  ierr = PetscContainerGetPointer(moabdata, (void**)&vmoab);CHKERRQ(ierr);
+
+  ierr = DMMoab_CreateVector(vmoab->mbiface,vmoab->pcomm,0,vmoab->tag_size,vmoab->ltog_tag,vmoab->tag_range,vmoab->serial,PETSC_TRUE,y);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoab_CreateTagName"
+/*  DMMoab_CreateTagName
+ *
+ *  Creates a unique tag name that will be shared across processes. If
+ *  pcomm is NULL, then this is a serial vector. A unique tag name
+ *  will be returned in tag_name in either case.
+ *
+ *  The tag names have the format _PETSC_VEC_N where N is some integer.
+ */
+PetscErrorCode DMMoab_CreateTagName(const moab::ParallelComm *pcomm,std::string& tag_name)
+{
+  moab::ErrorCode mberr;
+  PetscErrorCode  ierr;
+
+  PetscFunctionBegin;
+  const std::string PVEC_PREFIX      = "_PETSC_VEC_";
+  const PetscInt    PVEC_PREFIX_SIZE = PVEC_PREFIX.size();
+
+  // Check to see if there are any PETSc vectors defined:
+  const moab::Interface  *mbiface = pcomm->get_moab();
+  std::vector<moab::Tag> tags;
+  PetscInt               n = 0;
+  mberr = mbiface->tag_get_tags(tags);MBERRNM(mberr);
+  for(unsigned i = 0; i < tags.size(); i++) {
+    std::string s;
+    mberr = mbiface->tag_get_name(tags[i],s);MBERRNM(mberr);
+    if(s.find(PVEC_PREFIX) != std::string::npos){
+      // This tag represents a PETSc vector. Find the vector number:
+      PetscInt m;
+      std::istringstream(s.substr(PVEC_PREFIX_SIZE)) >> m;
+      if(m >= n) n = m+1;
+    }
+  }
+
+  // Make sure that n is consistent across all processes:
+  PetscInt global_n;
+  MPI_Comm comm = PETSC_COMM_SELF;
+  if(pcomm) comm = pcomm->comm();
+  ierr = MPI_Allreduce(&n,&global_n,1,MPI_INT,MPI_MAX,comm);CHKERRQ(ierr);
+
+  // Set the answer and return:
+  std::ostringstream ss;
+  ss << PVEC_PREFIX << global_n;
+  tag_name = ss.str();
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "DMMoab_VecUserDestroy"
+PetscErrorCode DMMoab_VecUserDestroy(void *user)
+{
+  Vec_MOAB        *vmoab;
+  PetscErrorCode  ierr;
+  moab::ErrorCode merr;
+
+  PetscFunctionBegin;
+  vmoab = (Vec_MOAB*)user;
+  vmoab->tag_range.~Range();
+  if(vmoab->new_tag) {
+    // Tag created via a call to VecDuplicate, delete the underlying tag in MOAB...
+    merr = vmoab->mbiface->tag_delete(vmoab->tag);MBERRNM(merr);
+  }
+
+  ierr = PetscFree(vmoab);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+}
+

File src/dm/impls/moab/examples/makefile

+
+ALL:
+
+LOCDIR	 = src/dm/impls/moab/examples/
+DIRS	 = tests
+
+include ${PETSC_DIR}/conf/variables
+include ${PETSC_DIR}/conf/rules
+include ${PETSC_DIR}/conf/test

File src/dm/impls/moab/examples/tests/ex1.cxx

+static char help[] = "Simple MOAB example\n\n";
+
+#include <petscdmmoab.h>
+#include <iostream>
+#include "moab/Interface.hpp"
+#include "moab/ScdInterface.hpp"
+#include "MBTagConventions.hpp"
+
+class AppCtx {
+public:
+  DM            dm;                /* REQUIRED in order to use SNES evaluation functions */
+  PetscLogEvent createMeshEvent;
+  /* Domain and mesh definition */
+  PetscInt dim;
+  char filename[PETSC_MAX_PATH_LEN];
+  char tagname[PETSC_MAX_PATH_LEN];
+  moab::Interface *iface;
+  moab::ParallelComm *pcomm;
+  moab::Tag tag;
+  moab::Tag ltog_tag;
+  moab::Range range;
+  PetscInt tagsize;
+
+  AppCtx()
+          : dm(NULL), dim(-1), iface(NULL), pcomm(NULL), tag(0), ltog_tag(0), tagsize(0)
+      {strcpy(filename, ""); strcpy(tagname, "");}
+
+};
+
+#undef __FUNCT__
+#define __FUNCT__ "ProcessOptions"
+PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  strcpy(options->filename, "");
+  strcpy(options->tagname, "");
+  options->dim = -1;
+
+  ierr = PetscOptionsBegin(comm, "", "MOAB example options", "DMMOAB");CHKERRQ(ierr);
+  ierr = PetscOptionsInt("-dim", "The topological mesh dimension", "ex1.c", options->dim, &options->dim, NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsString("-filename", "The file containing the mesh", "ex1.c", options->filename, options->filename, sizeof(options->filename), NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsString("-tagname", "The tag name from which to create a vector", "ex1.c", options->tagname, options->tagname, sizeof(options->tagname), NULL);CHKERRQ(ierr);
+  ierr = PetscOptionsEnd();
+
+  ierr = PetscLogEventRegister("CreateMesh",          DM_CLASSID,   &options->createMeshEvent);CHKERRQ(ierr);
+  PetscFunctionReturn(0);
+};
+
+#undef __FUNCT__
+#define __FUNCT__ "CreateMesh"
+PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
+{
+  PetscErrorCode ierr;
+
+  PetscFunctionBegin;
+  ierr = PetscLogEventBegin(user->createMeshEvent,0,0,0,0);CHKERRQ(ierr);
+  ierr = DMMoabCreateMoab(comm, user->iface, user->pcomm, user->ltog_tag, &user->range, dm);CHKERRQ(ierr);
+  std::cout << "Created DMMoab using DMMoabCreateDMAndInstance." << std::endl;
+  ierr = DMMoabGetInterface(*dm, &user->iface);CHKERRQ(ierr);
+
+    // load file and get entities of requested or max dimension
+  moab::ErrorCode merr;
+  if (strlen(user->filename) > 0) {
+    merr = user->iface->load_file(user->filename);MBERRNM(merr);
+    std::cout << "Read mesh from file " << user->filename << std::endl;
+  }
+  else {
+      // make a simple structured mesh
+    moab::ScdInterface *scdi;
+    merr = user->iface->query_interface(scdi);
+    moab::ScdBox *box;
+    merr = scdi->construct_box (moab::HomCoord(0,0,0), moab::HomCoord(2,2,2), NULL, 0, box);MBERRNM(merr);
+    user->dim = 3;
+    std::cout << "Created structured 2x2x2 mesh." << std::endl;
+  }
+  if (-1 == user->dim) {
+    moab::Range tmp_range;
+    merr = user->iface->get_entities_by_handle(0, tmp_range);MBERRNM(merr);
+    if (tmp_range.empty()) {
+      MBERRNM(moab::MB_FAILURE);
+    }
+    user->dim = user->iface->dimension_from_handle(*tmp_range.rbegin());
+  }
+  merr = user->iface->get_entities_by_dimension(0, user->dim, user->range);MBERRNM(merr);
+  ierr = DMMoabSetRange(*dm, user->range);CHKERRQ(ierr);
+
+    // get the requested tag if a name was input
+  if (strlen(user->tagname)) {
+    merr = user->iface->tag_get_handle(user->tagname, user->tag);MBERRNM(merr);
+    moab::DataType ttype;
+    merr = user->iface->tag_get_data_type(user->tag, ttype);MBERRNM(merr);
+    if (ttype != moab::MB_TYPE_DOUBLE) {
+      printf("Tag type must be double!.\n");
+      return 1;
+    }
+  }
+  else {
+      // make a new tag
+    merr = user->iface->tag_get_handle("petsc_tag", 1, moab::MB_TYPE_DOUBLE, user->tag, moab::MB_TAG_CREAT | moab::MB_TAG_DENSE); MBERRNM(merr);
+      // initialize new tag with gids
+    std::vector<double> tag_vals(user->range.size());
+    moab::Tag gid_tag;
+    merr = user->iface->tag_get_handle("GLOBAL_ID", gid_tag);MBERRNM(merr);
+    merr = user->iface->tag_get_data(gid_tag, user->range, tag_vals.data());MBERRNM(merr); // read them into ints
+    double *dval = tag_vals.data(); int *ival = reinterpret_cast<int*>(dval); // "stretch" them into doubles, from the end
+    for (int i = tag_vals.size()-1; i >= 0; i--) dval[i] = ival[i];
+    merr = user->iface->tag_set_data(user->tag, user->range, tag_vals.data());MBERRNM(merr); // write them into doubles
+  }
+  merr = user->iface->tag_get_length(user->tag, user->tagsize);MBERRNM(merr);
+
+    // create the dmmoab and initialize its data
+  ierr = PetscObjectSetName((PetscObject) *dm, "MOAB mesh");CHKERRQ(ierr);
+  ierr = PetscLogEventEnd(user->createMeshEvent,0,0,0,0);CHKERRQ(ierr);
+  user->dm = *dm;
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc, char **argv)
+{
+  AppCtx         user;                 /* user-defined work context */
+  PetscErrorCode ierr;
+  Vec            vec;
+
+  ierr = PetscInitialize(&argc, &argv, NULL, help);CHKERRQ(ierr);
+  ierr = ProcessOptions(PETSC_COMM_WORLD, &user);CHKERRQ(ierr);
+
+  ierr = CreateMesh(PETSC_COMM_WORLD, &user, &user.dm);CHKERRQ(ierr); /* create the MOAB dm and the mesh */
+  ierr = DMMoabCreateVector(user.dm, user.tag, 1, user.range, PETSC_TRUE, PETSC_FALSE,
+                              &vec);CHKERRQ(ierr); /* create a vec from user-input tag */
+  std::cout << "Created VecMoab from existing tag." << std::endl;
+  ierr = VecDestroy(&vec);CHKERRQ(ierr);
+  std::cout << "Destroyed VecMoab." << std::endl;
+  ierr = DMDestroy(&user.dm);CHKERRQ(ierr);
+  std::cout << "Destroyed DMMoab." << std::endl;
+  ierr = PetscFinalize();
+  return 0;
+}

File src/dm/impls/moab/examples/tests/makefile

+
+CFLAGS	        =
+FFLAGS	        =
+CPPFLAGS        =
+FPPFLAGS        =
+LOCDIR          = src/dm/impls/moab/examples/tests
+EXAMPLESC       = ex1.cxx
+EXAMPLESF       =
+MANSEC          = DM
+
+include ${PETSC_DIR}/conf/variables
+include ${PETSC_DIR}/conf/rules
+
+ex1: ex1.o  chkopts
+	-${CLINKER} -o ex1 ex1.o ${PETSC_DM_LIB}
+	${RM} -f ex1.o
+
+#--------------------------------------------------------------------------
+runex1:
+	-@${MPIEXEC} -n 1 ./ex1 > ex1_1.tmp 2>&1;\
+	   if (${DIFF} output/ex1_1.out ex1_1.tmp) then true ;  \
+	   else echo ${PWD} ; echo "Possible problem with with runex1, diffs above \n========================================="; fi ;\
+	   ${RM} -f ex1_1.tmp
+
+TESTEXAMPLES_MOAB = ex1.PETSc runex1 ex1.rm
+
+include ${PETSC_DIR}/conf/test

File src/dm/impls/moab/examples/tests/output/ex1_1.out

+Created DMMoab using DMMoabCreateDMAndInstance.
+Created structured 2x2x2 mesh.
+Created VecMoab from existing tag.
+Destroyed VecMoab.
+Destroyed DMMoab.

File src/dm/impls/moab/makefile

+#requirespackage  'PETSC_HAVE_MOAB'
+#requiresprecision double
+#requiresscalar    real
+
+ALL: lib
+
+CFLAGS   =
+FFLAGS   =
+SOURCEC  = dmmoab.cxx
+SOURCEF  =
+SOURCEH  =
+DIRS     = examples
+LIBBASE  = libpetscdm
+MANSEC   = DM
+LOCDIR   = src/dm/impls/moab/
+
+include ${PETSC_DIR}/conf/variables
+include ${PETSC_DIR}/conf/rules
+include ${PETSC_DIR}/conf/test
+

File src/dm/interface/dmregall.c

 PETSC_EXTERN PetscErrorCode DMCreate_Mesh(DM);
 PETSC_EXTERN PetscErrorCode DMCreate_Cartesian(DM);
 #endif
+PETSC_EXTERN PetscErrorCode  DMCreate_Moab(DM);
 
 #undef __FUNCT__
 #define __FUNCT__ "DMRegisterAll"
   ierr = DMRegisterDynamic(DMMESH,      path, "DMCreate_Mesh",      DMCreate_Mesh);CHKERRQ(ierr);
   ierr = DMRegisterDynamic(DMCARTESIAN, path, "DMCreate_Cartesian", DMCreate_Cartesian);CHKERRQ(ierr);
 #endif
+#if defined(PETSC_HAVE_MOAB)
+  ierr = DMRegisterDynamic(DMMOAB,      path, "DMCreate_Moab",      DMCreate_Moab);CHKERRQ(ierr);
+#endif
   PetscFunctionReturn(0);
 }
 

File src/ts/examples/tutorials/ex30.cxx

+static const char help[] = "Time-dependent Brusselator reaction-diffusion PDE in 1d. Demonstrates IMEX methods and uses MOAB.\n";
+/*
+   u_t - alpha u_xx = A + u^2 v - (B+1) u
+   v_t - alpha v_xx = B u - u^2 v
+   0 < x < 1;
+   A = 1, B = 3, alpha = 1/50
+
+   Initial conditions:
+   u(x,0) = 1 + sin(2 pi x)
+   v(x,0) = 3
+
+   Boundary conditions:
+   u(0,t) = u(1,t) = 1
+   v(0,t) = v(1,t) = 3
+*/
+
+// PETSc includes:
+#include <petscts.h>
+#include <petscdmmoab.h>
+
+// MOAB includes:
+#if defined (PETSC_HAVE_MOAB)
+#  include <moab/Core.hpp>
+#  include <moab/ReadUtilIface.hpp>
+#  include <MBTagConventions.hpp>
+#else
+#error You must have MOAB for this example. Reconfigure using --download-moab
+#endif
+
+typedef struct {
+  PetscScalar u,v;
+} Field;
+
+typedef struct _User *User;
+struct _User {
+  PetscReal A,B;                /* Reaction coefficients */
+  PetscReal alpha;              /* Diffusion coefficient */
+  PetscReal uleft,uright;       /* Dirichlet boundary conditions */
+  PetscReal vleft,vright;       /* Dirichlet boundary conditions */
+  PetscInt  npts;               /* Number of mesh points */
+
+  moab::ParallelComm *pcomm;
+  moab::Interface *mbint;
+  moab::Range *owned_vertexes;
+  moab::Range *owned_edges;
+  moab::Range *all_vertexes;
+  moab::Range *shared_vertexes;
+  moab::Tag    unknowns_tag;
+  PetscInt     unknowns_tag_size;
+  moab::Tag    id_tag;
+};
+
+static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
+static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
+static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat*,Mat*,MatStructure*,void*);
+
+PetscErrorCode create_app_data(_User& user);
+PetscErrorCode destroy_app_data(_User& user);
+PetscErrorCode create_matrix(_User &user, DM &dm, Mat *J);
+
+/****************
+ *              *
+ *     MAIN     *
+ *              *
+ ****************/
+#undef __FUNCT__
+#define __FUNCT__ "main"
+int main(int argc,char **argv)
+{
+  TS                ts;         /* nonlinear solver */
+  Vec               X;          /* solution, residual vectors */
+  Mat               J;          /* Jacobian matrix */
+  PetscInt          steps,maxsteps,mx;
+  PetscErrorCode    ierr;
+  PetscReal         ftime,dt;
+  _User             user;       /* user-defined work context */
+  TSConvergedReason reason;
+  DM                dm;
+  moab::ErrorCode   merr;
+
+  PetscInitialize(&argc,&argv,(char *)0,help);
+
+  // Fill in the user defined work context (creates mesh, solution field on mesh)
+  ierr = create_app_data(user);CHKERRQ(ierr);
+
+    // Create the DM to manage the mesh
+  ierr = DMMoabCreateMoab(user.pcomm->comm(),user.mbint,user.pcomm,user.id_tag,user.owned_vertexes,&dm);CHKERRQ(ierr);
+
+  // Create the solution vector:
+  ierr = DMMoabCreateVector(dm,user.unknowns_tag,user.unknowns_tag_size,*user.owned_vertexes,PETSC_FALSE,PETSC_FALSE,&X);CHKERRQ(ierr);
+    // Create the matrix
+  ierr = create_matrix(user, dm, &J);CHKERRQ(ierr);
+
+  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+     Create timestepping solver context
+     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+  ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
+  ierr = TSSetType(ts,TSARKIMEX);CHKERRQ(ierr);
+
+  ierr = TSSetRHSFunction(ts,PETSC_NULL,FormRHSFunction,&user);CHKERRQ(ierr);
+  ierr = TSSetIFunction(ts,PETSC_NULL,FormIFunction,&user);CHKERRQ(ierr);
+  ierr = TSSetIJacobian(ts,J,J,FormIJacobian,&user);CHKERRQ(ierr);
+
+  ftime = 10.0;
+  maxsteps = 10000;
+  ierr = TSSetDuration(ts,maxsteps,ftime);CHKERRQ(ierr);
+
+  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+     Set initial conditions
+   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+  ierr = TSSetSolution(ts,X);CHKERRQ(ierr);
+  ierr = VecGetSize(X,&mx);CHKERRQ(ierr);
+  dt   = 0.4 * user.alpha / PetscSqr(mx); /* Diffusive CFL */
+  ierr = TSSetInitialTimeStep(ts,0.0,dt);CHKERRQ(ierr);
+
+  // /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+  //    Set runtime options
+  //  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+  ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
+
+  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+     Solve nonlinear system
+     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+  ierr = TSSolve(ts,X);CHKERRQ(ierr);
+  ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
+  ierr = TSGetTimeStepNumber(ts,&steps);CHKERRQ(ierr);
+  ierr = TSGetConvergedReason(ts,&reason);CHKERRQ(ierr);
+  ierr = PetscPrintf(PETSC_COMM_WORLD,"%s at time %G after %D steps\n",TSConvergedReasons[reason],ftime,steps);CHKERRQ(ierr);
+
+
+  // Print the solution:
+  ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
+
+  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+     Write out the final mesh
+   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+  merr = user.mbint->write_file("ex30-final.h5m");MBERRNM(merr);
+
+  /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+     Free work space.
+   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
+
+  // Free all PETSc related resources:
+  ierr = MatDestroy(&J);CHKERRQ(ierr);
+  ierr = VecDestroy(&X);CHKERRQ(ierr);
+  ierr = TSDestroy(&ts);CHKERRQ(ierr);
+  ierr = DMDestroy(&dm);CHKERRQ(ierr);
+
+  // Free all MOAB related resources:
+  ierr = destroy_app_data(user);CHKERRQ(ierr);
+
+  ierr = PetscFinalize();
+  return 0;
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "FormIFunction"
+static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ptr)
+{
+  User            user = (User)ptr;
+  PetscInt        i;
+  Field           *x,*xdot,*f;
+  PetscReal       hx;
+  PetscErrorCode  ierr;
+  moab::ErrorCode merr;
+  int             *vertex_ids;
+
+  PetscFunctionBegin;
+  hx = 1.0/(PetscReal)(user->npts-1);
+
+  // Get connectivity information:
+  int verts_per_entity;
+  int count,num_edges;
+  moab::EntityHandle *connect;
+  moab::Range::iterator iter = user->owned_edges->begin();
+  merr = user->mbint->connect_iterate(iter,user->owned_edges->end(),connect,verts_per_entity,num_edges);MBERRNM(merr);
+
+  // Get tag data:
+  moab::Tag x_tag;
+  ierr = DMMoabGetVecTag(X,&x_tag);CHKERRQ(ierr);
+  merr = user->mbint->tag_iterate(x_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
+				  count,reinterpret_cast<void*&>(x));MBERRNM(merr);
+
+  moab::Tag xdot_tag;
+  ierr = DMMoabGetVecTag(Xdot,&xdot_tag);CHKERRQ(ierr);
+  merr = user->mbint->tag_iterate(xdot_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
+				  count,reinterpret_cast<void*&>(xdot));MBERRNM(merr);
+
+  moab::Tag f_tag;
+  ierr = DMMoabGetVecTag(F,&f_tag);CHKERRQ(ierr);
+  merr = user->mbint->tag_iterate(f_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
+				  count,reinterpret_cast<void*&>(f));MBERRNM(merr);
+
+  // Get the global IDs on all vertexes:
+  moab::Tag id_tag;
+  ierr = user->mbint->tag_get_handle(GLOBAL_ID_TAG_NAME,id_tag);CHKERRQ(ierr);
+  merr = user->mbint->tag_iterate(id_tag,user->all_vertexes->begin(),user->all_vertexes->end(),
+  				  count,reinterpret_cast<void*&>(vertex_ids));MBERRNM(merr);
+
+  // Exchange tags that are needed for assembly:
+  merr = user->pcomm->exchange_tags(x_tag,*user->all_vertexes);MBERRNM(merr);
+  merr = user->pcomm->exchange_tags(xdot_tag,*user->all_vertexes);MBERRNM(merr);
+
+  // Compute f:
+  const Field zero_field = {0.0,0.0};
+  std::fill(f,f+user->all_vertexes->size(),zero_field);
+
+  const moab::EntityHandle first_vertex = *user->all_vertexes->begin();
+
+  for (i = 0; i < num_edges; i++) {
+    const moab::EntityHandle idx_left  = connect[2*i]-first_vertex;
+    const moab::EntityHandle idx_right = connect[2*i+1]-first_vertex;
+
+    const int id_left  = vertex_ids[idx_left ];
+    const int id_right = vertex_ids[idx_right];
+
+    if (id_left == 0) {
+      // Apply left BC
+      f[idx_left].u += hx * (x[idx_left].u - user->uleft);
+      f[idx_left].v += hx * (x[idx_left].v - user->vleft);
+    } else {
+      f[idx_left].u += hx * xdot[idx_left].u - user->alpha*(-2*x[idx_left].u + x[idx_right].u)/hx;
+      f[idx_left].v += hx * xdot[idx_left].v - user->alpha*(-2*x[idx_left].v + x[idx_right].v)/hx;
+    }
+
+    if (id_right == user->npts-1) {
+      // Apply right BC
+      f[idx_right].u += hx * (x[idx_right].u - user->uright);
+      f[idx_right].v += hx * (x[idx_right].v - user->vright);
+    } else {
+      f[idx_right].u -= user->alpha*x[idx_left].u/hx;
+      f[idx_right].v -= user->alpha*x[idx_left].v/hx;
+    }
+  }
+
+  // Add tags on shared vertexes:
+  merr = user->pcomm->reduce_tags(f_tag,MPI_SUM,*user->shared_vertexes);MBERRNM(merr);
+
+  // Print vectors for debugging:
+  // ierr = PetscPrintf(PETSC_COMM_WORLD, "X\n");CHKERRQ(ierr);
+  // ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
+
+  // ierr = PetscPrintf(PETSC_COMM_WORLD, "Xdot\n");CHKERRQ(ierr);
+  // ierr = VecView(Xdot,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
+
+  // ierr = PetscPrintf(PETSC_COMM_WORLD, "F\n");CHKERRQ(ierr);
+  // ierr = VecView(F,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
+
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "FormRHSFunction"
+static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
+{
+  User            user = (User)ptr;
+  PetscReal       hx;
+  Field           *x,*f;
+  PetscErrorCode  ierr;
+  moab::Tag       id_tag;
+  moab::ErrorCode merr;
+
+  PetscFunctionBegin;
+  hx = 1.0/(PetscReal)(user->npts-1);
+
+  merr = user->mbint->tag_get_handle(GLOBAL_ID_TAG_NAME,id_tag);MBERRNM(merr);
+
+  /* Get pointers to vector data */
+  ierr = VecGetArray(X,reinterpret_cast<PetscScalar**>(&x));CHKERRQ(ierr);
+  ierr = VecGetArray(F,reinterpret_cast<PetscScalar**>(&f));CHKERRQ(ierr);
+
+  /* Compute function over the locally owned part of the grid */
+  const moab::EntityHandle first_vertex = *user->owned_vertexes->begin();
+  moab::Range::iterator iter;
+  for(iter = user->owned_vertexes->begin(); iter != user->owned_vertexes->end(); iter++) {
+    moab::EntityHandle i = *iter - first_vertex;
+    PetscScalar u = x[i].u,v = x[i].v;
+    f[i].u = hx*(user->A + u*u*v - (user->B+1)*u);
+    f[i].v = hx*(user->B*u - u*u*v);
+  }
+
+  /* Restore vectors */
+  ierr = VecRestoreArray(X,reinterpret_cast<PetscScalar**>(&x));CHKERRQ(ierr);
+  ierr = VecRestoreArray(F,reinterpret_cast<PetscScalar**>(&f));CHKERRQ(ierr);
+
+  // Print vectors for debugging:
+  /* ierr = PetscPrintf(PETSC_COMM_WORLD,"RHS");CHKERRQ(ierr); */
+  /* ierr = VecView(F,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
+
+  PetscFunctionReturn(0);
+}
+
+/* --------------------------------------------------------------------- */
+/*
+  IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
+*/
+#undef __FUNCT__
+#define __FUNCT__ "FormIJacobian"
+PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat *J,Mat *Jpre,MatStructure *str,void *ptr)
+{
+  User            user = (User)ptr;
+  PetscErrorCode  ierr;
+  PetscInt        i;
+  PetscReal       hx;
+  moab::Tag       id_tag;
+  moab::ErrorCode merr;
+
+  PetscFunctionBegin;
+  hx = 1.0/(PetscReal)(user->npts-1);
+
+  merr = user->mbint->tag_get_handle(GLOBAL_ID_TAG_NAME,id_tag);MBERRNM(merr);
+
+  /* Compute function over the locally owned part of the grid */
+  moab::Range::iterator iter;
+  for(iter = user->owned_vertexes->begin(); iter != user->owned_vertexes->end(); iter++) {
+    merr = user->mbint->tag_get_data(id_tag,&(*iter),1,&i);MBERRNM(merr);
+
+    if (i == 0 || i == user->npts-1) {
+      // Boundary conditions...
+      const PetscInt row = i,col = i;
+      const PetscScalar vals[2][2] = {{hx,0},{0,hx}};
+      ierr = MatSetValuesBlocked(*Jpre,1,&row,1,&col,&vals[0][0],INSERT_VALUES);CHKERRQ(ierr);
+    } else {
+      //
+      const PetscInt row = i,col[] = {i-1,i,i+1};
+      const PetscScalar dxxL = -user->alpha/hx,dxx0 = 2.*user->alpha/hx,dxxR = -user->alpha/hx;
+      const PetscScalar vals[2][3][2] = {{{dxxL,0},{a*hx+dxx0,0},{dxxR,0}},
+                                         {{0,dxxL},{0,a*hx+dxx0},{0,dxxR}}};
+      ierr = MatSetValuesBlocked(*Jpre,1,&row,3,col,&vals[0][0][0],INSERT_VALUES);CHKERRQ(ierr);
+    }
+  }
+
+  ierr = MatAssemblyBegin(*Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  ierr = MatAssemblyEnd(*Jpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  if (*J != *Jpre) {
+    ierr = MatAssemblyBegin(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+    ierr = MatAssemblyEnd(*J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
+  }
+
+  /* ierr = MatView(*J,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
+
+  PetscFunctionReturn(0);
+}
+
+
+/********************************
+ *                              *
+ *     initialize_moab_mesh     *
+ *                              *
+ ********************************/
+
+#undef __FUNCT__
+#define __FUNCT__ "initialize_moab_mesh"
+PetscErrorCode initialize_moab_mesh(moab::ParallelComm* pcomm,int npts,int nghost,moab::Tag &unknowns_tag,PetscInt &unknowns_tag_size,moab::Tag &id_tag)
+{
+  moab::ErrorCode merr;
+  PetscInt num_procs;
+  PetscInt rank;
+
+  PetscFunctionBegin;
+  MPI_Comm_size( pcomm->comm(),&num_procs );
+  MPI_Comm_rank( pcomm->comm(),&rank );
+
+  // Begin with some error checking:
+  if(npts < 2) {
+    SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"npts must be >= 2");
+  }
+
+  if(num_procs >= npts) {
+    SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_OUTOFRANGE,"Num procs must be < npts");
+  }
+
+  // No errors,proceed with building the mesh:
+  moab::Interface *mbint = pcomm->get_moab();
+  moab::ReadUtilIface* readMeshIface;
+  mbint->query_interface(readMeshIface);
+
+  // Determine which elements (cells) this process owns:
+  const PetscInt nele = npts-1;
+  PetscInt my_nele = nele / num_procs; // Number of elements owned by this process
+  const PetscInt extra = nele % num_procs;
+  PetscInt vstart = rank * my_nele;    // The starting element for this process
+  if(rank < extra) my_nele++;
+  vstart += std::min(extra, rank);
+  const PetscInt my_npts = my_nele + 1;
+
+  // Create the local portion of the mesh:
+
+  // Create vertexes and set the coodinate of each vertex:
+  moab::EntityHandle vertex_first;
+  std::vector<double*> vertex_coords;
+  const int sequence_size = (my_nele + 2) + 1;
+  merr = readMeshIface->get_node_coords(3,my_npts,1,vertex_first,vertex_coords,sequence_size);MBERRNM(merr);
+
+  const double xs = 0.0, xe = 1.0;
+  const double dx = (xe - xs) / nele;
+  for (int i = 0; i < my_npts; i++) {
+    vertex_coords[0][i] = (i+vstart)*dx;
+    vertex_coords[1][i] = vertex_coords[2][i] = 0.0;
+  }
+
+  moab::Range owned_vertexes;
+  merr = mbint->get_entities_by_type(0,moab::MBVERTEX,owned_vertexes);MBERRNM(merr);
+
+  // Create elements between mesh points. This is done so that VisIt
+  // will interpret the output as a mesh that can be plotted...
+  moab::EntityHandle edge_first;
+  moab::EntityHandle *connectivity = 0;
+  merr = readMeshIface->get_element_connect(my_nele,2,moab::MBEDGE,1,edge_first,connectivity);MBERRNM(merr);
+
+  for (int i = 0; i < my_nele; i+=1) {
+    connectivity[2*i]   = vertex_first + i;
+    connectivity[2*i+1] = vertex_first + (i+1);
+  }
+
+  merr = readMeshIface->update_adjacencies(edge_first,my_nele,2,connectivity);MBERRNM(merr);
+
+  // Set tags on all of the vertexes...
+
+  // Create tag handle to represent the unknowns, u and v and
+  // initialize them. We will create a single tag whose type is an
+  // array of two doubles and whose name is "unknowns"
+  Field default_value = {0.0,0.0};
+  unknowns_tag_size = sizeof(Field)/sizeof(PetscScalar);
+  merr = mbint->tag_get_handle("unknowns",unknowns_tag_size,moab::MB_TYPE_DOUBLE,unknowns_tag,
+                              moab::MB_TAG_DENSE | moab::MB_TAG_CREAT,&default_value);MBERRNM(merr);
+
+  // Apply the "unknowns" tag to the vertexes with some initial value...
+  std::vector<Field> tag_data(my_npts);
+  for (int i = 0; i < my_npts; i++) {
+    tag_data[i].u = 1+sin(2*PETSC_PI*vertex_coords[0][i]);
+    tag_data[i].v = 3;
+  }
+
+  merr = mbint->tag_set_data(unknowns_tag,owned_vertexes,tag_data.data());MBERRNM(merr);
+
+  // Get the global ID tag. The global ID tag is applied to each
+  // vertex. It acts as an global identifier which MOAB uses to
+  // assemble the individual pieces of the mesh:
+  merr = mbint->tag_get_handle(GLOBAL_ID_TAG_NAME,id_tag);MBERRNM(merr);
+
+  std::vector<int> global_ids(my_npts);
+  for (int i = 0; i < my_npts; i++) {
+    global_ids[i] = i+vstart;
+  }
+
+  merr = mbint->tag_set_data(id_tag,owned_vertexes,global_ids.data());MBERRNM(merr);
+
+  moab::Range owned_edges;
+  merr = mbint->get_entities_by_type(0,moab::MBEDGE,owned_edges);MBERRNM(merr);
+  merr = pcomm->resolve_shared_ents(0,owned_edges,1,0);MBERRNM(merr);
+
+  // Reassign global IDs on all entities.
+  merr = pcomm->assign_global_ids(0,1,0,false);MBERRNM(merr);
+  merr = pcomm->exchange_ghost_cells( 1,0,nghost,0,true);MBERRNM(merr);
+
+  // Everything is set up, now just do a tag exchange to update tags
+  // on all of the shared vertexes:
+  merr = pcomm->exchange_tags(id_tag,owned_vertexes);MBERRNM(merr);
+
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "create_app_data"
+PetscErrorCode create_app_data(_User& user)
+{
+  PetscErrorCode ierr;
+  moab::ErrorCode merr;
+
+  PetscFunctionBegin;
+
+  ierr = PetscOptionsBegin(PETSC_COMM_WORLD,PETSC_NULL,"Advection-reaction options",""); {
+    user.A      = 1;
+    user.B      = 3;
+    user.alpha  = 0.02;
+    user.uleft  = 1;
+    user.uright = 1;
+    user.vleft  = 3;
+    user.vright = 3;
+    user.npts   = 11;
+    ierr = PetscOptionsReal("-A","Reaction rate","",user.A,&user.A,PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscOptionsReal("-B","Reaction rate","",user.B,&user.B,PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscOptionsReal("-alpha","Diffusion coefficient","",user.alpha,&user.alpha,PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscOptionsReal("-uleft","Dirichlet boundary condition","",user.uleft,&user.uleft,PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscOptionsReal("-uright","Dirichlet boundary condition","",user.uright,&user.uright,PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscOptionsReal("-vleft","Dirichlet boundary condition","",user.vleft,&user.vleft,PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscOptionsReal("-vright","Dirichlet boundary condition","",user.vright,&user.vright,PETSC_NULL);CHKERRQ(ierr);
+    ierr = PetscOptionsInt("-npts","Number of mesh points","",user.npts,&user.npts,PETSC_NULL);CHKERRQ(ierr);
+  } ierr = PetscOptionsEnd();CHKERRQ(ierr);
+
+  user.mbint = new moab::Core;
+  user.pcomm = new moab::ParallelComm(user.mbint,PETSC_COMM_WORLD);
+
+  ierr = initialize_moab_mesh(user.pcomm,user.npts,1,user.unknowns_tag,user.unknowns_tag_size,user.id_tag); CHKERRQ(ierr); // creates moab mesh and field of unknowns on vertices
+
+  // Set the edge/vertex range once so we don't have to do it
+  // again. To do this we get all of the edges/vertexes then filter so
+  // we have only the owned entities in the ranges:
+  user.owned_edges = new moab::Range;
+  merr = user.mbint->get_entities_by_type(0,moab::MBEDGE,*user.owned_edges);MBERRNM(merr);
+
+  user.owned_vertexes = new moab::Range;
+  merr = user.mbint->get_entities_by_type(0,moab::MBVERTEX,*user.owned_vertexes);MBERRNM(merr);
+  user.all_vertexes = new moab::Range(*user.owned_vertexes);
+
+  user.shared_vertexes = new moab::Range;
+  merr = user.mbint->get_entities_by_type(0,moab::MBVERTEX,*user.shared_vertexes);MBERRNM(merr);
+
+  merr = user.pcomm->filter_pstatus(*user.owned_edges,PSTATUS_NOT_OWNED,PSTATUS_NOT);MBERRNM(merr);
+
+  merr = user.pcomm->filter_pstatus(*user.owned_vertexes,PSTATUS_NOT_OWNED,PSTATUS_NOT);MBERRNM(merr);
+
+  merr = user.pcomm->filter_pstatus(*user.shared_vertexes,PSTATUS_SHARED,PSTATUS_OR);MBERRNM(merr);
+
+  // Do some error checking...make sure that tag_data is in a single
+  // sequence:
+  int count;
+  void *data;
+  moab::Tag tag;
+  merr = user.mbint->tag_get_handle("unknowns",tag);MBERRNM(merr);
+  merr = user.mbint->tag_iterate(tag,user.all_vertexes->begin(),user.all_vertexes->end(),
+				 count,data);MBERRNM(merr);
+  if((unsigned int) count != user.all_vertexes->size()) {
+    SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_LIB,"Tag data not laid out contiguously %i %i",
+	     count,user.all_vertexes->size());
+  }
+
+  PetscFunctionReturn(0);
+}
+
+
+#undef __FUNCT__
+#define __FUNCT__ "create_matrix"
+PetscErrorCode create_matrix(_User &user, DM &dm, Mat *J)
+{
+  PetscErrorCode ierr;
+  moab::ErrorCode merr;
+  moab::Tag ltog_tag;
+  moab::Range range;
+
+  PetscFunctionBegin;
+  ierr = DMMoabGetLocalToGlobalTag(dm,&ltog_tag);CHKERRQ(ierr);
+  ierr = DMMoabGetRange(dm,&range);CHKERRQ(ierr);
+
+  // Create the matrix:
+  ierr = MatCreateBAIJ(user.pcomm->comm(),2,2*range.size(),2*range.size(),
+  		       PETSC_DECIDE,PETSC_DECIDE,3, NULL, 3, NULL, J);CHKERRQ(ierr);
+
+  // Set local to global numbering using the ltog_tag:
+  ISLocalToGlobalMapping ltog;
+  PetscInt               *gindices = new PetscInt[range.size()];
+  merr = user.pcomm->get_moab()->tag_get_data(ltog_tag, range, gindices);MBERRNM(merr);
+
+  ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, range.size(), gindices,PETSC_COPY_VALUES, &ltog);CHKERRQ(ierr);
+  ierr = MatSetLocalToGlobalMappingBlock(*J,ltog,ltog);CHKERRQ(ierr);
+
+  // Clean up:
+  ierr = ISLocalToGlobalMappingDestroy(&ltog);CHKERRQ(ierr);
+  delete [] gindices;
+
+  PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "destroy_app_data"
+PetscErrorCode destroy_app_data(_User& user)
+{
+  PetscFunctionBegin;
+
+  delete user.owned_edges;
+  delete user.owned_vertexes;
+  delete user.all_vertexes;
+  delete user.shared_vertexes;
+  delete user.pcomm;
+  delete user.mbint;
+
+  PetscFunctionReturn(0);
+}

File src/ts/examples/tutorials/makefile

 LOCDIR          = src/ts/examples/tutorials/
 EXAMPLESC       = ex1.c ex2.c ex3.c ex4.c ex5.c ex6.c ex7.c ex8.c \
                 ex9.c ex10.c ex12.c ex13.c ex14.c ex15.c ex16.c ex17.c \
-                ex19.c ex20.c ex21.c ex22.c ex23.c ex24.c ex25.c ex26.c ex28.c
+                ex19.c ex20.c ex21.c ex22.c ex23.c ex24.c ex25.c ex26.c ex28.c ex30.cxx
 EXAMPLESF       = ex1f.F ex2f.F ex22f.F # ex22f_mf.F90
 EXAMPLESFH      = ex2f.h
 MANSEC          = TS
 	-${CLINKER} -o ex29 ex29.o ${PETSC_TS_LIB}
 	${RM} ex29.o
 
+ex30: ex30.o chkopts
+	-${CLINKER} -o ex30 ex30.o ${PETSC_TS_LIB}
+	${RM} ex30.o
+
 #---------------------------------------------------------------------------------
 runex1:
 	-@${MPIEXEC} -n 1 ./ex1 -ksp_gmres_cgs_refinement_type refine_always -snes_type newtonls -ts_monitor_pseudo > ex1_1.tmp 2>&1;	  \
 	   ${DIFF} output/ex29.out ex29.tmp || echo  ${PWD} "\nPossible problem with ex29, diffs above \n=========================================";  \
 	   ${RM} -f ex29.tmp
 
+runex30:
+	-@${MPIEXEC} -n 1 ./ex30 > ex30_1.tmp 2>&1;	  \
+	   ${DIFF} output/ex30_1.out ex30_1.tmp || echo  ${PWD} "\nPossible problem with ex30_1, diffs above \n========================================="; \
+	   ${RM} -f ex30_1.tmp
 
 TESTEXAMPLES_C		  = ex1.PETSc runex1 runex1_2 ex1.rm ex2.PETSc runex2 ex2.rm ex3.PETSc runex3 runex3_2 ex3.rm \
                             ex4.PETSc runex4 runex4_2 runex4_3 runex4_4 ex4.rm ex5.PETSc ex5.rm \
 TESTEXAMPLES_C_X_MPIUNI =
 TESTEXAMPLES_13		  = ex2.PETSc ex2.rm ex3.PETSc ex3.rm ex4.PETSc ex4.rm \
                             ex5.PETSc ex5.rm
-TESTEXAMPLES_EXODUSII = ex11.PETSc runex11 ex11.rm
+TESTEXAMPLES_EXODUSII     = ex11.PETSc runex11 ex11.rm
+TESTEXAMPLES_MOAB         = ex30.PETSc runex30 ex30.rm
 
 include ${PETSC_DIR}/conf/test

File src/ts/examples/tutorials/output/ex30_1.out

+CONVERGED_TIME at time 10.1839 after 51 steps
+Vector Object: 1 MPI processes
+  type: mpi
+Process [0]
+1
+3
+1.61629
+2.10951
+1.88038
+1.69013
+1.88146
+1.65113
+1.66999
+1.91526
+1.24455
+2.45621
+0.791781
+3.02322
+0.54165
+3.34375
+0.508752
+3.42344
+0.668931
+3.30342
+1
+3

File src/ts/examples/tutorials/output/ex30_2.out

+CONVERGED_TIME at time 10.1839 after 51 steps
+Vector Object: 2 MPI processes
+  type: mpi
+Process [0]
+1
+3
+1.61629
+2.10951
+1.88038
+1.69013
+1.88146
+1.65113
+1.66999
+1.91526
+1.24455
+2.45621
+Process [1]
+0.791781
+3.02322
+0.54165
+3.34375
+0.508752
+3.42344
+0.668931
+3.30342
+1
+3