Commits

BarryFSmith committed 7e25b11 Merge with conflicts

Merge branch 'knepley/remove-sieve' into next

I guessed about the ViennaCL stuff

Conflicts:
include/petscsnes.h
src/dm/impls/mesh/examples/tests/output/ex10_0.out
src/dm/impls/mesh/examples/tests/output/ex10_1.out
src/dm/impls/mesh/examples/tests/output/ex11f90_0.out
src/dm/impls/mesh/examples/tests/output/ex11f90_1.out
src/dm/impls/mesh/examples/tests/output/ex3_0.out
src/dm/impls/mesh/examples/tests/output/ex3_1.out
src/dm/impls/mesh/examples/tests/output/ex3_2.out
src/dm/impls/mesh/examples/tests/output/ex3_3.out
src/dm/impls/mesh/examples/tests/output/ex8_0.out
src/dm/impls/mesh/examples/tests/output/ex9f90_0.out
src/dm/impls/mesh/examples/tests/output/ex9f90_1.out
src/dm/impls/mesh/mesh.c
src/dm/interface/dlregisdmdm.c
src/vec/vec/interface/vecregall.c

Comments (0)

Files changed (288)

bin/processSummary.py

-#!/usr/bin/env python
-import user
-import script
-
-class LogSummarizer(script.Script):
-  transforms = [
-    (', boost::tuples::null_type', ''),
-    ('ALE::pair<ALE::Point, int>', 'Rank'),
-    ('ALE::Sieve<ALE::Point, int, int>', 'Topology'),
-    ('ALE::CoSifter<Topology, ALE::Point, ALE::Point, int>', 'Bundle'),
-    ('ALE::CoSifter<Topology, ALE::Point, ALE::Point, double>', 'Field'),
-    ('ALE::CoSifter<Topology, Rank, ALE::Point, int>', 'RankBundle'),
-    ('ALE::CoSifter<Topology, Rank, ALE::Point, double>', 'RankField'),
-    ('std::less<ALE::Point>', 'lessPoint'),
-    ('std::less<Rank >', 'lessRank'),
-    ('Bundle::trueFunc<ALE::Point>', 'trueFunc'),
-    ('RankField::trueFunc<ALE::Point>', 'trueFunc'),
-    ('Field::trueFunc<ALE::Point>', 'trueFunc'),
-    ('boost::multi_index::composite_key_compare<lessPoint, lessPoint, lessPoint>', 'stdCompare'),
-    ('boost::multi_index::composite_key_compare<lessPoint, std::less<int>, lessPoint>', 'stdCompareInt'),
-    ('boost::multi_index::composite_key_compare<lessPoint, trueFunc, lessPoint>', 'orderCompare'),
-    ('boost::multi_index::composite_key_compare<lessPoint, trueFunc, lessRank>', 'orderRankCompare'),
-    ('ALE::array<ALE::Point>', 'PointArray'),
-    ('ALE::set<ALE::Point>', 'PointSet'),
-    ('std::vector<ALE::Point, std::allocator<ALE::Point> >', 'PointVec'),
-    ('std::set<int, std::less<int>, std::allocator<int> >', 'IntSet'),
-    ('ALE::SifterDef::Arrow<ALE::Point, ALE::Point, ALE::Point>', 'Arrow'),
-    ('ALE::SifterDef::Arrow<ALE::Point, ALE::Point, int>', 'IntArrow'),
-    ('ALE::SifterDef::Arrow<ALE::Point, Rank, ALE::Point>', 'RankArrow'),
-    ('ALE::SifterDef::Arrow<int, ALE::Point, ALE::pair<ALE::Point, ALE::pair<int, int> > >', 'OverlapArrow'),
-    ('ALE::SifterDef::Arrow<ALE::Point, int, ALE::pair<ALE::Point, ALE::pair<int, int> > >', 'FlipOverlapArrow'),
-    ('ALE::SifterDef::Arrow<int, ALE::pair<int, ALE::Point>, ALE::pair<ALE::Point, ALE::pair<int, int> > >', 'BiOverlapArrow'),
-    ('ALE::SifterDef::Arrow<ALE::pair<int, ALE::Point>, int, ALE::pair<ALE::Point, ALE::pair<int, int> > >', 'FlipBiOverlapArrow'),
-    ('ALE::SifterDef::Rec<ALE::Point>', 'Point'),
-    ('ALE::SifterDef::Rec<Rank >', 'RankPoint'),
-    ('ALE::SifterDef::Rec<int>', 'IntPoint'),
-    ('ALE::SifterDef::Rec<ALE::pair<int, ALE::Point> >', 'IPPairPoint'),
-    ('ALE::SieveDef::Rec<ALE::Point, int>', 'Point&Int'),
-    ('ALE::SieveDef::Rec<int, ALE::Point>', 'Int&Point'),
-    ('ALE::SifterDef::RecContainerTraits<ALE::Point, Point >', 'SifterRecTraits'),
-    ('ALE::SifterDef::RecContainerTraits<int, IntPoint >', 'OverlapRecTraits'),
-    ('ALE::SieveDef::RecContainerTraits<ALE::Point, Point&Int >', ' SieveRecTraits'),
-    ('ALE::SifterDef::RecContainer<ALE::Point, Point >', 'SifterRecCon'),
-    ('ALE::SifterDef::RecContainer<Rank, RankPoint >', 'RankRecCon'),
-    ('ALE::SifterDef::RecContainer<int, IntPoint >', 'OverlapRecCon'),
-    ('ALE::SifterDef::RecContainer<ALE::pair<int, ALE::Point>, IPPairPoint >', 'BiOverlapRecCon'),
-    ('ALE::SieveDef::RecContainer<ALE::Point, Point&Int >', 'SieveRecCon'),
-    ('SifterRecTraits::PointSequence', 'SifPointSeq'),
-    ('OverlapRecTraits::PointSequence', 'OvPointSeq'),
-    ('SieveRecTraits::PointSequence', 'SivPointSeq'),
-    ('SieveRecTraits::TwoValueSequence< SieveRecTraits::heightMarkerTag, int>', 'HeightSeq'),
-    ('SieveRecTraits::TwoValueSequence< SieveRecTraits::depthMarkerTag, int>', 'DepthSeq'),
-    ('ALE::RightSequenceDuplicator<ALE::ConeArraySequence<Arrow > >', 'SeqDup'),
-    ('ALE::RightSequenceDuplicator<ALE::ConeArraySequence<IntArrow > >', 'IntSeqDup'),
-    ('ALE::Sifter<ALE::Point, ALE::Point, ALE::Point, orderCompare, SifterRecCon, SifterRecCon >', 'Order'),
-    ('ALE::Sifter<ALE::Point, ALE::Point, ALE::Point, stdCompare, SifterRecCon, SifterRecCon >', 'ReOrder'),
-    ('ALE::Sifter<ALE::Point, Rank, ALE::Point, orderRankCompare, SifterRecCon, RankRecCon >', 'RankOrder'),
-    ('ALE::ASifter<int, ALE::Point, ALE::pair<ALE::Point, ALE::pair<int, int> >, (ALE::SifterDef::ColorMultiplicity)1, boost::multi_index::composite_key_compare<std::less<int>, std::less<ALE::pair<ALE::Point, ALE::pair<int, int> > >, lessPoint>, OverlapRecCon, SifterRecCon >', 'Overlap'),
-    ('ALE::ASifter<ALE::Point, int, ALE::pair<ALE::Point, ALE::pair<int, int> >, (ALE::SifterDef::ColorMultiplicity)1, boost::multi_index::composite_key_compare<lessPoint, std::less<ALE::pair<ALE::Point, ALE::pair<int, int> > >, std::less<int>>, SifterRecCon, OverlapRecCon >', 'FlipOverlap'),
-    ('ALE::ASifter<int, ALE::pair<int, ALE::Point>, ALE::pair<ALE::Point, ALE::pair<int, int> >, (ALE::SifterDef::ColorMultiplicity)1, boost::multi_index::composite_key_compare<std::less<int>, std::less<ALE::pair<ALE::Point, ALE::pair<int, int> > >, std::less<ALE::pair<int, ALE::Point> >>, OverlapRecCon, BiOverlapRecCon >', 'BiOverlap'),
-    ('ALE::ASifter<ALE::pair<int, ALE::Point>, int, ALE::pair<ALE::Point, ALE::pair<int, int> >, (ALE::SifterDef::ColorMultiplicity)1, boost::multi_index::composite_key_compare<std::less<ALE::pair<int, ALE::Point> >, std::less<ALE::pair<ALE::Point, ALE::pair<int, int> > >, std::less<int>>, BiOverlapRecCon, OverlapRecCon >', 'FlipBiOverlap'),
-    ('ALE::ASifter<ALE::Point, ALE::Point, int, (ALE::SifterDef::ColorMultiplicity)1, stdCompareInt, SieveRecCon, SieveRecCon >::traits::coneSequence', 'ConeSeq'),
-    ('ALE::ASifter<ALE::Point, ALE::Point, int, (ALE::SifterDef::ColorMultiplicity)1, stdCompareInt, SieveRecCon, SieveRecCon >::traits::supportSequence', 'SuppSeq'),
-    ('ALE::ASifter<ALE::Point, ALE::Point, ALE::Point, (ALE::SifterDef::ColorMultiplicity)1, stdCompare, SifterRecCon, SifterRecCon >::traits::coneSequence', 'OrderFusionCone'),
-    ('ALE::ASifter<ALE::Point, ALE::Point, ALE::Point, (ALE::SifterDef::ColorMultiplicity)1, stdCompare, SifterRecCon, SifterRecCon >::traits::supportSequence', 'OrderFusionSupp'),
-    ('ALE::ASifter<ALE::Point, ALE::Point, ALE::Point, (ALE::SifterDef::ColorMultiplicity)1, orderCompare, SifterRecCon, SifterRecCon >::traits::coneSequence', 'OrderCone'),
-    ('ALE::ASifter<ALE::Point, ALE::Point, ALE::Point, (ALE::SifterDef::ColorMultiplicity)1, orderCompare, SifterRecCon, SifterRecCon >::traits::supportSequence', 'OrderSupp'),
-    ('BiOverlap::traits::coneSequence', 'BiOverlapCone'),
-    ('BiOverlap::traits::supportSequence', 'BiOverlapSupp'),
-    ('FlipBiOverlap::traits::coneSequence', 'FlipBiOverlapCone'),
-    ('FlipBiOverlap::traits::supportSequence', 'FlipBiOverlapSupp'),
-    ('ALE::Flip<Order >', 'FlipOrder'),
-    ('ALE::Flip<ReOrder >', 'FlipReOrder'),
-    ('boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<Arrow > > > >', 'bArrow'),
-    ('boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<IntArrow > > > >', 'bIntArrow'),
-    ('boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<RankArrow > > > >', 'bRankArrow'),
-    ('boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<OverlapArrow > > > >', 'bOverlapArrow'),
-    ('boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<FlipOverlapArrow > > > >', 'bFlipOverlapArrow'),
-    ('boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<BiOverlapArrow > > > >', 'bBiOverlapArrow'),
-    ('boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<FlipBiOverlapArrow > > > >', 'bFlipBiOverlapArrow'),
-    ('boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<Point > >', 'bPoint'),
-    ('boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<Point&Int > > > > >', 'bIntPoint'),
-    ('boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<RankPoint > >', 'bRankPoint'),
-    ('boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<OverlapPoint > > > >', 'bOverlapPoint'),
-    ('boost::multi_index::detail::copy_map_entry<bArrow >', 'bArrowCopy'),
-    ('boost::multi_index::detail::copy_map_entry<bPoint >', 'bPointCopy'),
-    ('std::_Rb_tree_node<ALE::Point>', 'tPoint')
-    ]
-
-  removes = [
-    'Registered event',
-    'boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<ALE::SifterDef::Arrow<ALE::Point, ALE::Point, int> > > > >:',
-    'boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<ALE::SieveDef::Rec<ALE::Point, int> > > > > >:',
-    'boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<ALE::SifterDef::Arrow<ALE::Point, ALE::Point, ALE::Point> > > > >:',
-    'boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<ALE::SifterDef::Rec<ALE::Point> > >:',
-    'boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<ALE::SifterDef::Arrow<ALE::Point, ALE::pair<ALE::Point, int>, ALE::Point> > > > >:',
-    'boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<ALE::SifterDef::Rec<ALE::pair<ALE::Point, int> > > >:',
-    'boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<ALE::SifterDef::Rec<int> > >',
-    'boost::multi_index::detail::ordered_index_node<boost::multi_index::detail::index_node_base<ALE::SifterDef::Rec<ALE::pair<int, ALE::Point> > > >',
-    'std::_Rb_tree_node<ALE::Point>:',
-    'Obj:std::set<std::string',
-    'std::set<std::string'
-    ]
-
-  def setupHelp(self, help):
-    help = script.Script.setupHelp(self, help)
-    #help.addArgument('LogSummarizer', '-', nargs.ArgBool(None, 0, 'Print this help message', isTemporary = 1), ignoreDuplicates = 1)
-    return help
-
-  def filter(self, line):
-    for r in LogSummarizer.removes:
-      if line.startswith(r):
-        return ''
-    return line
-
-  def transform(self, line):
-    '''We could try to pick out just the event lines'''
-    for longName, shortName in LogSummarizer.transforms:
-      line = line.replace(longName, shortName)
-    return line
-
-  def process(self, inF, outF):
-    for line in inF.readlines():
-      line = self.filter(line)
-      line = self.transform(line)
-      outF.write(line)
-    return
-
-  def run(self):
-    import sys
-
-    self.setup()
-    self.process(sys.stdin, sys.stdout)
-    return
-
-if __name__ == '__main__':
-  LogSummarizer().run()

config/PETSc/packages/Sieve.py

-import PETSc.package
-
-class Configure(PETSc.package.NewPackage):
-  def __init__(self, framework):
-    import os
-    PETSc.package.NewPackage.__init__(self, framework)
-    self.include         = [os.path.abspath(os.path.join('include', 'sieve'))]
-    self.libdir          = ''
-    self.archIndependent = 1
-    self.cxx             = 1
-    self.complex         = 1
-    self.worksonWindows  = 1
-    self.double          = 0
-    return
-
-  def setupDependencies(self, framework):
-    PETSc.package.NewPackage.setupDependencies(self, framework)
-    self.boost    = self.framework.require('config.packages.boost',self)
-    self.mpi      = self.framework.require('config.packages.MPI',self)
-    #self.triangle = self.framework.require('PETSc.packages.Triangle',self)
-    #self.tetgen   = self.framework.require('PETSc.packages.TetGen',self)
-    self.deps = [self.boost, self.mpi]
-    return
-
-  def getSearchDirectories(self):
-    return [self.petscdir.dir]
-
-  def setupHelp(self, help):
-    PETSc.package.NewPackage.setupHelp(self, help)
-    import nargs
-    help.addArgument('Sieve', '-with-opt-sieve=<bool>', nargs.ArgBool(None, 1, 'Use IMesh which are optimized for interval point sets'))
-    help.addArgument('Sieve', '-with-sieve-memory-logging=<bool>', nargs.ArgBool(None, 0, 'Turn on memory logging for Sieve objects'))
-    return
-
-  def Install(self):
-    import sys
-    sieveDir = self.getDir()
-    self.framework.actions.addArgument('Sieve', 'Install', 'Installed Sieve into '+sieveDir)
-    return sieveDir
-
-  def consistencyChecks(self):
-    PETSc.package.NewPackage.consistencyChecks(self)
-    if self.framework.argDB['with-'+self.package]:
-      if not self.languages.clanguage == 'Cxx':
-        raise RuntimeError('Sieve requires C++. Suggest using --with-clanguage=cxx')
-      if not self.boost.found:
-        raise RuntimeError('Sieve requires boost, and configure could not locate it. Suggest using --download-boost=1')
-      if 'with-opt-sieve' in self.argDB and self.argDB['with-opt-sieve']:
-        self.addDefine('OPT_SIEVE', 1)
-        self.addDefine('MESH_TYPE', 'ALE::IMesh<PetscInt, PetscScalar>')
-      else:
-        self.addDefine('MESH_TYPE', 'ALE::Mesh<PetscInt, PetscScalar>')
-      if 'with-sieve-memory-logging' in self.argDB and self.argDB['with-sieve-memory-logging']:
-        self.framework.addDefine('ALE_MEM_LOGGING', 1)
-    return

config/examples/arch-linux-cxx-sieve-gcov.py

-#!/usr/bin/env python
-
-configure_options = [
-  '--with-cc=gcc',
-  '--with-cxx=g++',
-  '--download-mpich=1',
-  '--with-sieve=1',
-  '--download-boost=1',
-  '--with-clanguage=cxx',
-  #'--with-gcov=1',
-  ]
-
-if __name__ == '__main__':
-  import sys,os
-  sys.path.insert(0,os.path.abspath('config'))
-  import configure
-  configure.petsc_configure(configure_options)

config/examples/arch-linux-cxx-sieve.py

-#!/usr/bin/env python
-
-configure_options = [
-  '--with-cc=gcc',
-  '--with-cxx=g++',
-  '--download-mpich=1',
-  '--with-sieve=1',
-  '--download-boost=1',
-  '--with-clanguage=cxx'
-  ]
-
-if __name__ == '__main__':
-  import sys,os
-  sys.path.insert(0,os.path.abspath('config'))
-  import configure
-  configure.petsc_configure(configure_options)

include/finclude/ftn-custom/petscdmmesh.h90

-
-#if !defined(PETSC_USE_FORTRAN_MODULES)
-#include "finclude/ftn-custom/petscdmmeshdef.h90"
-#endif
-
-#include "finclude/ftn-custom/petscdmhide.h90"
-
-#if defined(PETSC_USE_FORTRAN_DATATYPES) && !defined(SECTIONREAL_HIDE)
-#define SECTIONREAL_HIDE type(SectionReal)
-#define SECTIONINT_HIDE type(SectionInt)
-#elif !defined(SECTIONREAL_HIDE)
-#define SECTIONREAL_HIDE SectionReal
-#define SECTIONINT_HIDE SectionInt
-#endif
-
-#if defined(PETSC_USE_FORTRAN_DATATYPES) && !defined(USE_DMMESH_HIDE)
-#  define USE_DMMESH_HIDE use petscdmmeshdef
-#elif !defined(USE_DMMESH_HIDE)
-#  define USE_DMMESH_HIDE
-#endif
-
-      Interface
-        Subroutine DMMeshGetCoordinatesF90(m,array,ierr)
-          USE_DMMESH_HIDE
-          PetscReal, pointer :: array(:,:)
-          PetscErrorCode ierr
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine DMMeshRestoreCoordinatesF90(m,array,ierr)
-          USE_DMMESH_HIDE
-          PetscReal, pointer :: array(:,:)
-          PetscErrorCode ierr
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine DMMeshGetElementsF90(m,array,ierr)
-          USE_DMMESH_HIDE
-          PetscInt, pointer :: array(:,:)
-          PetscErrorCode ierr
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine DMMeshRestoreElementsF90(m,array,ierr)
-          USE_DMMESH_HIDE
-          PetscInt, pointer :: array(:,:)
-          PetscErrorCode ierr
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine DMMeshGetConeF90(m,p,array,ierr)
-          USE_DMMESH_HIDE
-          PetscInt, pointer :: array(:)
-          PetscInt p
-          PetscErrorCode ierr
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine DMMeshRestoreConeF90(m,p,array,ierr)
-          USE_DMMESH_HIDE
-          PetscInt, pointer :: array(:)
-          PetscInt p
-          PetscErrorCode ierr
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine SectionGetArrayF90(m,name,array,ierr)
-          USE_DMMESH_HIDE
-          CHARACTER*80 name
-          PetscReal, pointer :: array(:,:)
-          PetscErrorCode ierr
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine SectionGetArray1DF90(m,name,array,ierr)
-          USE_DMMESH_HIDE
-          CHARACTER*80 name
-          PetscReal, pointer :: array(:)
-          PetscErrorCode ierr
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine BCSectionGetArrayF90(m,name,array,ierr)
-          USE_DMMESH_HIDE
-          CHARACTER*80 name
-          PetscInt, pointer :: array(:,:)
-          PetscErrorCode ierr
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine BCSectionGetArray1DF90(m,name,array,ierr)
-          USE_DMMESH_HIDE
-          CHARACTER*80 name
-          PetscInt, pointer :: array(:)
-          PetscErrorCode ierr
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine BCSectionRealGetArrayF90(m,name,array,ierr)
-          USE_DMMESH_HIDE
-          CHARACTER*80 name
-          PetscReal, pointer :: array(:,:)
-          PetscErrorCode ierr
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine BCFUNCGetArrayF90(m,array,ierr)
-          USE_DMMESH_HIDE
-          PetscReal, pointer :: array(:,:)
-          PetscErrorCode ierr
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine DMMeshGetLabelIds(m,name,array,ierr)
-          USE_DMMESH_HIDE
-          PetscErrorCode ierr
-          PetscInt, pointer :: array(:)
-          CHARACTER*80 name
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine DMMeshGetStratum(m,name,value,array,ierr)
-          USE_DMMESH_HIDE
-          PetscErrorCode ierr
-          PetscInt, pointer :: array(:)
-          PetscInt value
-          CHARACTER*80 name
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine SectionRealRestrict(section,e,values,ierr)
-          USE_DMMESH_HIDE
-          PetscErrorCode ierr
-          PetscScalar, pointer :: values(:)
-          PetscInt e
-          SECTIONREAL_HIDE section
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine SectionIntRestrict(section,e,values,ierr)
-          USE_DMMESH_HIDE
-          PetscErrorCode ierr
-          PetscInt, pointer :: values(:)
-          PetscInt e
-          SECTIONINT_HIDE section
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine SectionRealRestore(section,e,values,ierr)
-          USE_DMMESH_HIDE
-          PetscErrorCode ierr
-          PetscScalar, pointer :: values(:)
-          PetscInt e
-          SECTIONREAL_HIDE section
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine SectionIntRestore(section,e,values,ierr)
-          USE_DMMESH_HIDE
-          PetscErrorCode ierr
-          PetscInt, pointer :: values(:)
-          PetscInt e
-          SECTIONINT_HIDE section
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine SectionRealRestrictClosure(section,m,e,n,values,ierr)
-          USE_DMMESH_HIDE
-          PetscErrorCode ierr
-          PetscScalar, pointer :: values(:)
-          PetscInt n, e
-          SECTIONREAL_HIDE section
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine SectionIntRestrictClosure(section,m,e,n,values,ierr)
-          USE_DMMESH_HIDE
-          PetscErrorCode ierr
-          PetscInt, pointer :: values(:)
-          PetscInt n, e
-          SECTIONINT_HIDE section
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine SectionRealUpdate(section,e,values,mode,ierr)
-          USE_DMMESH_HIDE
-          PetscErrorCode ierr
-          PetscScalar, pointer :: values(:)
-          PetscInt e
-          InsertMode mode
-          SECTIONREAL_HIDE section
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine SectionIntUpdate(section,e,values,mode,ierr)
-          USE_DMMESH_HIDE
-          PetscErrorCode ierr
-          PetscInt, pointer :: values(:)
-          PetscInt e
-          InsertMode mode
-          SECTIONINT_HIDE section
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine SectionRealUpdateClosure(section,m,e,values,mode,    &
-     &ierr)
-          USE_DMMESH_HIDE
-          PetscErrorCode ierr
-          PetscScalar, pointer :: values(:)
-          PetscInt e
-          InsertMode mode
-          SECTIONREAL_HIDE section
-          DM_HIDE m
-        End Subroutine
-      End Interface
-
-      Interface
-        Subroutine SectionIntUpdateClosure(section,m,e,values,mode,ierr)
-          USE_DMMESH_HIDE
-          PetscErrorCode ierr
-          PetscInt, pointer :: values(:)
-          PetscInt e
-          InsertMode mode
-          SECTIONINT_HIDE section
-          DM_HIDE m
-        End Subroutine
-      End Interface

include/finclude/ftn-custom/petscdmmeshdef.h90

-#if defined(PETSC_USE_FORTRAN_DATATYPES)
-      type SectionReal
-        PetscFortranAddr:: v
-      end type SectionReal
-      type SectionInt
-        PetscFortranAddr:: v
-      end type SectionInt
-#endif

include/finclude/makefile

 FFLAGS    =
 SOURCEC   =
 SOURCEF   =
-SOURCEH   = petsc.h petscsys.h petsclog.h petscvec.h petscsnes.h petscdm.h petscdmda.h petscdmmesh.h petscdraw.h petscmat.h \
+SOURCEH   = petsc.h petscsys.h petsclog.h petscvec.h petscsnes.h petscdm.h petscdmda.h petscdraw.h petscmat.h \
 	    petscksp.h petscpc.h petscviewer.h petscis.h petscao.h \
 	    petscts.h petscis.h90 petscvec.h90 petscmat.h90 petscdm.h90 petscdmda.h90 petscdmadda.h90 petscdmcomposite.h90 \
-            petscdmmesh.h90 petscdmredundant.h90 \
-            petscdef.h  petscsysdef.h petsclogdef.h petscvecdef.h petscsnesdef.h petscdmdef.h petscdmdadef.h petscdmmeshdef.h \
+            petscdmredundant.h90 \
+            petscdef.h  petscsysdef.h petsclogdef.h petscvecdef.h petscsnesdef.h petscdmdef.h petscdmdadef.h \
             petscdrawdef.h petscmatdef.h petsckspdef.h petscpcdef.h petscviewerdef.h petscisdef.h petscaodef.h \
 	    petsctsdef.h
 LIBBASE   = libpetscvec

include/finclude/petsc.h

 #include "finclude/petscdmplex.h"
 !#include "finclude/petscdmadda.h"
 !#include "finclude/petscdmcomposite.h"
-#include "finclude/petscdmmesh.h"
 !#include "finclude/petscdmsliced.h"
 #include "finclude/petscts.h"

include/finclude/petsc.h90

 #include "finclude/petscdmadda.h90"
 #include "finclude/petscdmcomposite.h90"
 #include "finclude/petscdmredundant.h90"
-#include "finclude/petscdmmesh.h90"
 #include "finclude/petscpc.h90"
 #include "finclude/petscksp.h90"
 #include "finclude/petscsnes.h90"

include/finclude/petscdef.h

 #include "finclude/petscmatdef.h"
 #include "finclude/petscdmdef.h"
 #include "finclude/petscdmdadef.h"
-#include "finclude/petscdmmeshdef.h"
+#include "finclude/petscdmplexdef.h"
 #include "finclude/petscpcdef.h"
 #include "finclude/petsckspdef.h"
 #include "finclude/petscsnesdef.h"

include/finclude/petscdmmesh.h

-!
-!  Include file for Fortran use of the Mesh package in PETSc
-!
-#include "finclude/petscdmmeshdef.h"
-

include/finclude/petscdmmesh.h90

-!
-!
-!  Additional DMMesh include file for use of PETSc with Fortran 90/HPF
-!
-#include "finclude/ftn-custom/petscdmmesh.h90"
-#if defined(PETSC_USE_FORTRAN_INTERFACES)
-      interface
-#include "finclude/ftn-auto/petscdmmesh.h90"
-      end interface
-#endif

include/finclude/petscdmmeshdef.h

-!
-!
-!  Include file for Fortran use of the Mesh package in PETSc
-!
-#if !defined (__PETSCDMMESHDEF_H)
-#define __PETSCDMMESHDEF_H
-
-#include "finclude/petscdmdef.h"
-
-#if !defined(PETSC_USE_FORTRAN_DATATYPES)
-#define SectionReal PetscFortranAddr
-#define SectionInt  PetscFortranAddr
-#endif
-
-#endif

include/petsc-private/meshimpl.h

-#if !defined(_MESHIMPL_H)
-#define _MESHIMPL_H
-
-#include <petscmat.h>    /*I      "petscmat.h"    I*/
-#include <petscdmmesh.h> /*I      "petscdmmesh.h"    I*/
-#include "petsc-private/dmimpl.h"
-
-typedef struct Sieve_Label *SieveLabel;
-struct Sieve_Label {
-  char      *name;           /* Label name */
-  PetscInt   numStrata;      /* Number of integer values */
-  PetscInt  *stratumValues;  /* Value of each stratum */
-  PetscInt  *stratumOffsets; /* Offset of each stratum */
-  PetscInt  *stratumSizes;   /* Size of each stratum */
-  PetscInt  *points;         /* Points for each stratum, sorted after setup */
-  SieveLabel next;           /* Linked list */
-};
-
-typedef struct {
-  ALE::Obj<PETSC_MESH_TYPE> m;
-
-  PetscSection         defaultSection;
-  VecScatter           globalScatter;
-  DMMeshLocalFunction1 lf;
-  DMMeshLocalJacobian1 lj;
-
-  /*-------- NEW_MESH_IMPL -------------*/
-  PetscBool            useNewImpl;
-  PetscInt             dim; /* Topological mesh dimension */
-  PetscSF              sf;  /* SF for parallel point overlap */
-
-  /*   Sieve */
-  PetscSection         coneSection;    /* Layout of cones (inedges for DAG) */
-  PetscInt             maxConeSize;    /* Cached for fast lookup */
-  PetscInt            *cones;          /* Cone for each point */
-  PetscInt            *coneOrientations; /* TODO */
-  PetscSection         supportSection; /* Layout of cones (inedges for DAG) */
-  PetscInt             maxSupportSize; /* Cached for fast lookup */
-  PetscInt            *supports;       /* Cone for each point */
-  PetscSection         coordSection;   /* Layout for coordinates */
-  Vec                  coordinates;    /* Coordinate values */
-
-  PetscInt            *meetTmpA,    *meetTmpB;    /* Work space for meet operation */
-  PetscInt            *joinTmpA,    *joinTmpB;    /* Work space for join operation */
-  PetscInt            *closureTmpA, *closureTmpB; /* Work space for closure operation */
-
-  /* Labels */
-  SieveLabel           labels;         /* Linked list of labels */
-} DM_Mesh;
-
-typedef struct {
-  ALE::Obj<ALE::CartesianMesh> m;
-} DM_Cartesian;
-
-typedef struct _SectionRealOps *SectionRealOps;
-struct _SectionRealOps {
-  PetscErrorCode (*view)(SectionReal,PetscViewer);
-  PetscErrorCode (*restrictClosure)(SectionReal,int,PetscScalar**);
-  PetscErrorCode (*update)(SectionReal,int,const PetscScalar*,InsertMode);
-};
-
-struct _p_SectionReal {
-  PETSCHEADER(struct _SectionRealOps);
-  ALE::Obj<PETSC_MESH_TYPE::real_section_type> s;
-  ALE::Obj<PETSC_MESH_TYPE> b;
-};
-
-PETSC_EXTERN PetscClassId SECTIONREAL_CLASSID;
-PETSC_EXTERN PetscLogEvent SectionReal_View;
-
-typedef struct _SectionIntOps *SectionIntOps;
-struct _SectionIntOps {
-  PetscErrorCode (*view)(SectionInt,PetscViewer);
-  PetscErrorCode (*restrictClosure)(SectionInt,int,PetscInt**);
-  PetscErrorCode (*update)(SectionInt,int,const PetscInt*,InsertMode);
-};
-
-struct _p_SectionInt {
-  PETSCHEADER(struct _SectionIntOps);
-  ALE::Obj<PETSC_MESH_TYPE::int_section_type> s;
-  ALE::Obj<PETSC_MESH_TYPE> b;
-};
-
-PETSC_EXTERN PetscClassId SECTIONINT_CLASSID;
-PETSC_EXTERN PetscLogEvent SectionInt_View;
-
-#endif /* _MESHIMPL_H */
 #include <petscdmda.h>
 #include <petscdmadda.h>
 #include <petscdmcomposite.h>
-#include <petscdmmesh.h>
 #include <petscdmpatch.h>
 #include <petscdmplex.h>
 #include <petscdmredundant.h>

include/petscdm.h

 #define DMCOMPOSITE "composite"
 #define DMSLICED    "sliced"
 #define DMSHELL     "shell"
-#define DMMESH      "mesh"
 #define DMPLEX      "plex"
 #define DMCARTESIAN "cartesian"
 #define DMREDUNDANT "redundant"

include/petscdmmesh.h

-/*
-  DMMesh, for easy parallelism of simple unstructured distributed mesh problems.
-*/
-#if !defined(__PETSCDMMESH_H)
-#define __PETSCDMMESH_H
-
-#include <petscsftypes.h>
-#include <petscdm.h>
-
-#if defined(PETSC_HAVE_SIEVE) && defined(__cplusplus)
-
-#include <sieve/Mesh.hh>
-#include <sieve/CartesianSieve.hh>
-#include <sieve/Distribution.hh>
-#include <sieve/Generator.hh>
-
-/*S
-  DMMESH - DM object that encapsulates an unstructured mesh. This uses the Sieve package to represent the mesh.
-
-  Level: intermediate
-
-  Concepts: grids, grid refinement
-
-.seealso:  DM, DMMeshCreate()
-S*/
-/*S
-  DMCARTESIAN - DM object that encapsulates a structured, Cartesian mesh in any dimension. This is currently a
-  replacement for DMDA in dimensions greater than 3.
-
-  Level: intermediate
-
-  Concepts: grids, grid refinement
-
-.seealso:  DM, DMCartesianCreate(), DMDA
-S*/
-PETSC_EXTERN PetscLogEvent DMMesh_View, DMMesh_GetGlobalScatter, DMMesh_restrictVector, DMMesh_assembleVector, DMMesh_assembleVectorComplete, DMMesh_assembleMatrix, DMMesh_updateOperator;
-
-PETSC_EXTERN PetscErrorCode DMMeshCreate(MPI_Comm, DM*);
-PETSC_EXTERN PetscErrorCode DMMeshCreateMeshFromAdjacency(MPI_Comm, PetscInt, PetscInt, PetscInt, PetscInt [], PetscInt, PetscInt, const PetscReal [], PetscBool, DM*);
-PETSC_EXTERN PetscErrorCode DMMeshClone(DM, DM *);
-PETSC_EXTERN PetscErrorCode DMMeshGetMesh(DM, ALE::Obj<PETSC_MESH_TYPE>&);
-PETSC_EXTERN PetscErrorCode DMMeshSetMesh(DM, const ALE::Obj<PETSC_MESH_TYPE>&);
-PETSC_EXTERN PetscErrorCode DMMeshGetGlobalScatter(DM, VecScatter *);
-PETSC_EXTERN PetscErrorCode DMMeshFinalize();
-
-/* New Sieve Mesh interface */
-PETSC_EXTERN PetscErrorCode DMMeshGetDimension(DM, PetscInt *);
-PETSC_EXTERN PetscErrorCode DMMeshSetDimension(DM, PetscInt);
-PETSC_EXTERN PetscErrorCode DMMeshGetChart(DM, PetscInt *, PetscInt *);
-PETSC_EXTERN PetscErrorCode DMMeshSetChart(DM, PetscInt, PetscInt);
-PETSC_EXTERN PetscErrorCode DMMeshGetConeSize(DM, PetscInt, PetscInt *);
-PETSC_EXTERN PetscErrorCode DMMeshSetConeSize(DM, PetscInt, PetscInt);
-PETSC_EXTERN PetscErrorCode DMMeshGetCone(DM, PetscInt, const PetscInt *[]);
-PETSC_EXTERN PetscErrorCode DMMeshSetCone(DM, PetscInt, const PetscInt[]);
-PETSC_EXTERN PetscErrorCode DMMeshGetSupportSize(DM, PetscInt, PetscInt *);
-PETSC_EXTERN PetscErrorCode DMMeshGetSupport(DM, PetscInt, const PetscInt *[]);
-PETSC_EXTERN PetscErrorCode DMMeshGetConeSection(DM, PetscSection *);
-PETSC_EXTERN PetscErrorCode DMMeshGetCones(DM, PetscInt *[]);
-PETSC_EXTERN PetscErrorCode DMMeshGetMaxSizes(DM, PetscInt *, PetscInt *);
-PETSC_EXTERN PetscErrorCode DMMeshSetUp(DM);
-PETSC_EXTERN PetscErrorCode DMMeshSymmetrize(DM);
-PETSC_EXTERN PetscErrorCode DMMeshStratify(DM);
-
-PETSC_EXTERN PetscErrorCode DMMeshGetLabelValue(DM, const char[], PetscInt, PetscInt *);
-PETSC_EXTERN PetscErrorCode DMMeshSetLabelValue(DM, const char[], PetscInt, PetscInt);
-PETSC_EXTERN PetscErrorCode DMMeshGetLabelSize(DM, const char[], PetscInt *);
-PETSC_EXTERN PetscErrorCode DMMeshGetLabelIdIS(DM, const char[], IS *);
-PETSC_EXTERN PetscErrorCode DMMeshGetStratumSize(DM, const char [], PetscInt, PetscInt *);
-PETSC_EXTERN PetscErrorCode DMMeshGetStratumIS(DM, const char [], PetscInt, IS *);
-
-PETSC_EXTERN PetscErrorCode DMMeshJoinPoints(DM, PetscInt, const PetscInt [], PetscInt *);
-PETSC_EXTERN PetscErrorCode DMMeshMeetPoints(DM, PetscInt, const PetscInt [], PetscInt *, const PetscInt **);
-PETSC_EXTERN PetscErrorCode DMMeshGetTransitiveClosure(DM, PetscInt, PetscBool, PetscInt *, const PetscInt *[]);
-
-PETSC_EXTERN PetscErrorCode DMMeshCreatePartition(DM, PetscSection *, IS *, PetscInt);
-PETSC_EXTERN PetscErrorCode DMMeshCreatePartitionClosure(DM, PetscSection, IS, PetscSection *, IS *);
-
-/* Old Sieve Mesh interface */
-PETSC_EXTERN PetscErrorCode DMMeshDistribute(DM, const char[], DM*);
-PETSC_EXTERN PetscErrorCode DMMeshDistributeByFace(DM, const char[], DM*);
-PETSC_EXTERN PetscErrorCode DMMeshGenerate(DM, PetscBool , DM *);
-PETSC_EXTERN PetscErrorCode DMMeshRefine(DM, double, PetscBool , DM*);
-PETSC_EXTERN PetscErrorCode DMMeshLoad(PetscViewer, DM);
-PETSC_EXTERN PetscErrorCode DMMeshGetMaximumDegree(DM, PetscInt *);
-
-PETSC_EXTERN PetscErrorCode DMCartesianCreate(MPI_Comm, DM *);
-PETSC_EXTERN PetscErrorCode DMMeshCartesianGetMesh(DM, ALE::Obj<ALE::CartesianMesh>&);
-PETSC_EXTERN PetscErrorCode DMMeshCartesianSetMesh(DM, const ALE::Obj<ALE::CartesianMesh>&);
-
-PETSC_EXTERN PetscErrorCode DMMeshGetCoordinates(DM, PetscBool , PetscInt *, PetscInt *, PetscReal *[]);
-PETSC_EXTERN PetscErrorCode DMMeshGetElements(DM, PetscBool , PetscInt *, PetscInt *, PetscInt *[]);
-
-PETSC_EXTERN PetscErrorCode DMMeshCreateBoxMesh(MPI_Comm, PetscInt, PetscBool, DM *);
-PETSC_EXTERN PetscErrorCode DMMeshMarkBoundaryCells(DM, const char [], PetscInt, PetscInt);
-PETSC_EXTERN PetscErrorCode DMMeshGetDepthStratum(DM, PetscInt, PetscInt *, PetscInt *);
-PETSC_EXTERN PetscErrorCode DMMeshGetHeightStratum(DM, PetscInt, PetscInt *, PetscInt *);
-PETSC_EXTERN PetscErrorCode DMMeshCreateSection(DM, PetscInt, PetscInt, PetscInt [], PetscInt [], PetscInt, PetscInt [], IS [], PetscSection *);
-PETSC_EXTERN PetscErrorCode DMMeshConvertSection(const ALE::Obj<PETSC_MESH_TYPE>&, const ALE::Obj<PETSC_MESH_TYPE::real_section_type>&, PetscSection *);
-PETSC_EXTERN PetscErrorCode DMMeshGetSection(DM, const char [], PetscSection *);
-PETSC_EXTERN PetscErrorCode DMMeshSetSection(DM, const char [], PetscSection);
-PETSC_EXTERN PetscErrorCode DMMeshGetDefaultSection(DM, PetscSection *);
-PETSC_EXTERN PetscErrorCode DMMeshSetDefaultSection(DM, PetscSection);
-PETSC_EXTERN PetscErrorCode DMMeshGetCoordinateSection(DM, PetscSection *);
-PETSC_EXTERN PetscErrorCode DMMeshSetCoordinateSection(DM, PetscSection);
-PETSC_EXTERN PetscErrorCode DMMeshCreateConeSection(DM, PetscSection *);
-PETSC_EXTERN PetscErrorCode DMMeshGetCoordinateVec(DM, Vec *);
-PETSC_EXTERN PetscErrorCode DMMeshComputeCellGeometry(DM, PetscInt, PetscReal *, PetscReal *, PetscReal *, PetscReal *);
-PETSC_EXTERN PetscErrorCode DMMeshVecSetClosure(DM, Vec, PetscInt, const PetscScalar [], InsertMode);
-PETSC_EXTERN PetscErrorCode DMMeshVecGetClosure(DM, Vec, PetscInt, const PetscScalar **);
-PETSC_EXTERN PetscErrorCode DMMeshVecRestoreClosure(DM, Vec, PetscInt, const PetscScalar **);
-PETSC_EXTERN PetscErrorCode DMMeshMatSetClosure(DM, Mat, PetscInt, PetscScalar [], InsertMode);
-
-PETSC_EXTERN PetscErrorCode MatSetValuesTopology(Mat, DM, PetscInt, const PetscInt [], DM, PetscInt, const PetscInt [], const PetscScalar [], InsertMode);
-PETSC_EXTERN PetscErrorCode DMMeshRestrictVector(Vec, Vec, InsertMode);
-PETSC_EXTERN PetscErrorCode DMMeshAssembleVectorComplete(Vec, Vec, InsertMode);
-PETSC_EXTERN PetscErrorCode DMMeshAssembleVector(Vec, PetscInt, PetscScalar [], InsertMode);
-PETSC_EXTERN PetscErrorCode DMMeshUpdateOperator(Mat, const ALE::Obj<PETSC_MESH_TYPE>&, const ALE::Obj<PETSC_MESH_TYPE::real_section_type>&, const ALE::Obj<PETSC_MESH_TYPE::order_type>&, const PETSC_MESH_TYPE::point_type&, PetscScalar [], InsertMode);
-PETSC_EXTERN PetscErrorCode DMMeshUpdateOperatorGeneral(Mat, const ALE::Obj<PETSC_MESH_TYPE>&, const ALE::Obj<PETSC_MESH_TYPE::real_section_type>&, const ALE::Obj<PETSC_MESH_TYPE::order_type>&, const PETSC_MESH_TYPE::point_type&, const ALE::Obj<PETSC_MESH_TYPE>&, const ALE::Obj<PETSC_MESH_TYPE::real_section_type>&, const ALE::Obj<PETSC_MESH_TYPE::order_type>&, const PETSC_MESH_TYPE::point_type&, PetscScalar [], InsertMode);
-
-/*S
-  SectionReal - Abstract PETSc object that manages distributed field data over a topology (Sieve).
-
-  Level: beginner
-
-  Concepts: distributed mesh, field
-
-.seealso:  SectionRealCreate(), SectionRealDestroy(), Mesh, DMMeshCreate()
-S*/
-typedef struct _p_SectionReal* SectionReal;
-
-/* Logging support */
-PETSC_EXTERN PetscClassId SECTIONREAL_CLASSID;
-
-PETSC_EXTERN PetscErrorCode SectionRealCreate(MPI_Comm,SectionReal*);
-PETSC_EXTERN PetscErrorCode SectionRealDestroy(SectionReal*);
-PETSC_EXTERN PetscErrorCode SectionRealView(SectionReal,PetscViewer);
-PETSC_EXTERN PetscErrorCode SectionRealDuplicate(SectionReal,SectionReal*);
-
-PETSC_EXTERN PetscErrorCode SectionRealGetSection(SectionReal,ALE::Obj<PETSC_MESH_TYPE::real_section_type>&);
-PETSC_EXTERN PetscErrorCode SectionRealSetSection(SectionReal,const ALE::Obj<PETSC_MESH_TYPE::real_section_type>&);
-PETSC_EXTERN PetscErrorCode SectionRealGetBundle(SectionReal,ALE::Obj<PETSC_MESH_TYPE>&);
-PETSC_EXTERN PetscErrorCode SectionRealSetBundle(SectionReal,const ALE::Obj<PETSC_MESH_TYPE>&);
-
-PETSC_EXTERN PetscErrorCode SectionRealDistribute(SectionReal, DM, SectionReal *);
-PETSC_EXTERN PetscErrorCode SectionRealRestrict(SectionReal, PetscInt, PetscScalar *[]);
-PETSC_EXTERN PetscErrorCode SectionRealUpdate(SectionReal, PetscInt, const PetscScalar [], InsertMode);
-PETSC_EXTERN PetscErrorCode SectionRealZero(SectionReal);
-PETSC_EXTERN PetscErrorCode SectionRealCreateLocalVector(SectionReal, Vec*);
-PETSC_EXTERN PetscErrorCode SectionRealAddSpace(SectionReal);
-PETSC_EXTERN PetscErrorCode SectionRealGetFibration(SectionReal, const PetscInt, SectionReal *);
-PETSC_EXTERN PetscErrorCode SectionRealToVecDM(SectionReal, DM, ScatterMode, Vec);
-PETSC_EXTERN PetscErrorCode SectionRealToVec(SectionReal, VecScatter, ScatterMode, Vec);
-PETSC_EXTERN PetscErrorCode SectionRealNorm(SectionReal, DM, NormType, PetscReal *);
-PETSC_EXTERN PetscErrorCode SectionRealAXPY(SectionReal, DM, PetscScalar, SectionReal);
-PETSC_EXTERN PetscErrorCode SectionRealComplete(SectionReal);
-PETSC_EXTERN PetscErrorCode SectionRealSet(SectionReal, PetscReal);
-PETSC_EXTERN PetscErrorCode SectionRealGetFiberDimension(SectionReal, PetscInt, PetscInt*);
-PETSC_EXTERN PetscErrorCode SectionRealSetFiberDimension(SectionReal, PetscInt, const PetscInt);
-PETSC_EXTERN PetscErrorCode SectionRealSetFiberDimensionField(SectionReal, PetscInt, const PetscInt, const PetscInt);
-PETSC_EXTERN PetscErrorCode SectionRealGetSize(SectionReal, PetscInt *);
-PETSC_EXTERN PetscErrorCode SectionRealAllocate(SectionReal);
-PETSC_EXTERN PetscErrorCode SectionRealClear(SectionReal);
-
-PETSC_EXTERN PetscErrorCode SectionRealRestrictClosureWithArray(SectionReal, DM, PetscInt, PetscInt, PetscScalar []);
-PETSC_EXTERN PetscErrorCode SectionRealRestrictClosure(SectionReal, DM, PetscInt, const PetscScalar *[]);
-PETSC_EXTERN PetscErrorCode SectionRealUpdateClosure(SectionReal, DM, PetscInt, PetscScalar [], InsertMode);
-
-PETSC_EXTERN PetscErrorCode DMMeshHasSectionReal(DM, const char [], PetscBool  *);
-PETSC_EXTERN PetscErrorCode DMMeshGetSectionReal(DM, const char [], SectionReal *);
-PETSC_EXTERN PetscErrorCode DMMeshSetSectionReal(DM, const char [], SectionReal);
-PETSC_EXTERN PetscErrorCode DMMeshCreateVector(DM, SectionReal, Vec *);
-PETSC_EXTERN PetscErrorCode DMMeshCreateMatrix(DM, SectionReal, MatType, Mat *);
-PETSC_EXTERN PetscErrorCode DMMeshCreateGlobalScatter(DM, SectionReal, VecScatter *);
-PETSC_EXTERN PetscErrorCode DMMeshAssembleVectorDM(Vec, DM, SectionReal, PetscInt, PetscScalar [], InsertMode);
-PETSC_EXTERN PetscErrorCode DMMeshAssembleMatrixDM(Mat, DM, SectionReal, PetscInt, PetscScalar [], InsertMode);
-PETSC_EXTERN PetscErrorCode DMMeshSetupSection(DM, SectionReal);
-
-PETSC_EXTERN_TYPEDEF typedef PetscErrorCode (*DMMeshLocalFunction1)(DM, Vec, Vec, void*);
-PETSC_EXTERN_TYPEDEF typedef PetscErrorCode (*DMMeshLocalJacobian1)(DM, Vec, Mat, void*);
-
-PETSC_EXTERN PetscErrorCode DMMeshCreateGlobalRealVector(DM, SectionReal, Vec *);
-PETSC_EXTERN PetscErrorCode DMMeshGetGlobalScatter(DM,VecScatter *);
-PETSC_EXTERN PetscErrorCode DMMeshCreateGlobalScatter(DM,SectionReal,VecScatter *);
-PETSC_EXTERN PetscErrorCode DMMeshGetLocalFunction(DM, PetscErrorCode (**)(DM, Vec, Vec, void*));
-PETSC_EXTERN PetscErrorCode DMMeshSetLocalFunction(DM, PetscErrorCode (*)(DM, Vec, Vec, void*));
-PETSC_EXTERN PetscErrorCode DMMeshGetLocalJacobian(DM, PetscErrorCode (**)(DM, Vec, Mat, void*));
-PETSC_EXTERN PetscErrorCode DMMeshSetLocalJacobian(DM, PetscErrorCode (*)(DM, Vec, Mat, void*));
-
-/*S
-  SectionInt - Abstract PETSc object that manages distributed field data over a topology (Sieve).
-
-  Level: beginner
-
-  Concepts: distributed mesh, field
-
-.seealso:  SectionIntCreate(), SectionIntDestroy(), DM, DMMeshCreate()
-S*/
-typedef struct _p_SectionInt* SectionInt;
-
-/* Logging support */
-PETSC_EXTERN PetscClassId SECTIONINT_CLASSID;
-
-PETSC_EXTERN PetscErrorCode SectionIntCreate(MPI_Comm,SectionInt*);
-PETSC_EXTERN PetscErrorCode SectionIntDestroy(SectionInt*);
-PETSC_EXTERN PetscErrorCode SectionIntView(SectionInt,PetscViewer);
-
-PETSC_EXTERN PetscErrorCode SectionIntGetSection(SectionInt,ALE::Obj<PETSC_MESH_TYPE::int_section_type>&);
-PETSC_EXTERN PetscErrorCode SectionIntSetSection(SectionInt,const ALE::Obj<PETSC_MESH_TYPE::int_section_type>&);
-PETSC_EXTERN PetscErrorCode SectionIntGetBundle(SectionInt,ALE::Obj<PETSC_MESH_TYPE>&);
-PETSC_EXTERN PetscErrorCode SectionIntSetBundle(SectionInt,const ALE::Obj<PETSC_MESH_TYPE>&);
-
-PETSC_EXTERN PetscErrorCode SectionIntDistribute(SectionInt, DM, SectionInt *);
-PETSC_EXTERN PetscErrorCode SectionIntRestrict(SectionInt, PetscInt, PetscInt *[]);
-PETSC_EXTERN PetscErrorCode SectionIntUpdate(SectionInt, PetscInt, const PetscInt [], InsertMode);
-PETSC_EXTERN PetscErrorCode SectionIntZero(SectionInt);
-PETSC_EXTERN PetscErrorCode SectionIntComplete(SectionInt);
-PETSC_EXTERN PetscErrorCode SectionIntGetFiberDimension(SectionInt, PetscInt, PetscInt*);
-PETSC_EXTERN PetscErrorCode SectionIntSetFiberDimension(SectionInt, PetscInt, const PetscInt);
-PETSC_EXTERN PetscErrorCode SectionIntSetFiberDimensionField(SectionInt, PetscInt, const PetscInt, const PetscInt);
-PETSC_EXTERN PetscErrorCode SectionIntGetSize(SectionInt, PetscInt *);
-PETSC_EXTERN PetscErrorCode SectionIntAllocate(SectionInt);
-PETSC_EXTERN PetscErrorCode SectionIntClear(SectionInt);
-
-PETSC_EXTERN PetscErrorCode SectionIntAddSpace(SectionInt);
-PETSC_EXTERN PetscErrorCode SectionIntGetFibration(SectionInt, const PetscInt, SectionInt *);
-PETSC_EXTERN PetscErrorCode SectionIntSet(SectionInt, PetscInt);
-
-PETSC_EXTERN PetscErrorCode SectionIntRestrictClosureWithArray(SectionInt, DM, PetscInt, PetscInt, PetscInt []);
-PETSC_EXTERN PetscErrorCode SectionIntUpdateClosure(SectionInt, DM, PetscInt, PetscInt [], InsertMode);
-
-PETSC_EXTERN PetscErrorCode DMMeshHasSectionInt(DM, const char [], PetscBool  *);
-PETSC_EXTERN PetscErrorCode DMMeshGetSectionInt(DM, const char [], SectionInt *);
-PETSC_EXTERN PetscErrorCode DMMeshSetSectionInt(DM, SectionInt);
-
-/* Misc Mesh functions*/
-PETSC_EXTERN PetscErrorCode DMMeshSetMaxDof(DM, PetscInt);
-PETSC_EXTERN PetscErrorCode SectionGetArray(DM, const char [], PetscInt *, PetscInt *, PetscScalar *[]);
-
-/* Helper functions for simple distributions */
-PETSC_EXTERN PetscErrorCode DMMeshGetVertexSectionReal(DM, const char[], PetscInt, SectionReal *);
-PETSC_EXTERN PetscErrorCode DMMeshGetVertexSectionInt(DM, const char[], PetscInt, SectionInt *);
-PETSC_EXTERN PetscErrorCode DMMeshGetCellSectionReal(DM, const char[], PetscInt, SectionReal *);
-PETSC_EXTERN PetscErrorCode DMMeshGetCellSectionInt(DM, const char[], PetscInt, SectionInt *);
-PETSC_EXTERN PetscErrorCode DMMeshCreateSectionRealIS(DM,IS,const char [],PetscInt,SectionReal *);
-
-/* Scatter for simple distributions */
-PETSC_EXTERN PetscErrorCode DMMeshCreateScatterToZeroVertex(DM,VecScatter *);
-PETSC_EXTERN PetscErrorCode DMMeshCreateScatterToZeroVertexSet(DM,IS,IS,VecScatter *);
-PETSC_EXTERN PetscErrorCode DMMeshCreateScatterToZeroCell(DM,VecScatter *);
-PETSC_EXTERN PetscErrorCode DMMeshCreateScatterToZeroCellSet(DM,IS,IS,VecScatter *);
-
-/* Support for various mesh formats */
-PETSC_EXTERN PetscErrorCode DMMeshCreateExodus(MPI_Comm, const char [], DM *);
-PETSC_EXTERN PetscErrorCode DMMeshCreateExodusNG(MPI_Comm, PetscInt, DM *);
-PETSC_EXTERN PetscErrorCode DMMeshExodusGetInfo(DM, PetscInt *, PetscInt *, PetscInt *, PetscInt *, PetscInt *);
-PETSC_EXTERN PetscErrorCode DMMeshViewExodusSplit(DM,PetscInt);
-
-PETSC_EXTERN PetscErrorCode VecViewExodusVertex(DM,Vec,MPI_Comm,PetscInt,PetscInt,PetscInt);
-PETSC_EXTERN PetscErrorCode VecLoadExodusVertex(DM,Vec,MPI_Comm,PetscInt,PetscInt,PetscInt);
-PETSC_EXTERN PetscErrorCode VecViewExodusVertexSet(DM,Vec,PetscInt,MPI_Comm,PetscInt,PetscInt,PetscInt);
-PETSC_EXTERN PetscErrorCode VecLoadExodusVertexSet(DM,Vec,PetscInt,MPI_Comm,PetscInt,PetscInt,PetscInt);
-PETSC_EXTERN PetscErrorCode VecViewExodusCell(DM,Vec,MPI_Comm,PetscInt,PetscInt,PetscInt);
-PETSC_EXTERN PetscErrorCode VecLoadExodusCell(DM,Vec,MPI_Comm,PetscInt,PetscInt,PetscInt);
-PETSC_EXTERN PetscErrorCode VecViewExodusCellSet(DM,Vec,PetscInt,MPI_Comm,PetscInt,PetscInt,PetscInt);
-PETSC_EXTERN PetscErrorCode VecLoadExodusCellSet(DM,Vec,PetscInt,MPI_Comm,PetscInt,PetscInt,PetscInt);
-
-PETSC_EXTERN PetscErrorCode DMMeshCreatePCICE(MPI_Comm, const int, const char [], const char [], PetscBool , const char [], DM *);
-
-PETSC_EXTERN PetscErrorCode DMWriteVTKHeader(PetscViewer);
-PETSC_EXTERN PetscErrorCode DMWriteVTKVertices(DM, PetscViewer);
-PETSC_EXTERN PetscErrorCode DMWriteVTKElements(DM, PetscViewer);
-PETSC_EXTERN PetscErrorCode DMWritePCICEVertices(DM, PetscViewer);
-PETSC_EXTERN PetscErrorCode DMWritePCICEElements(DM, PetscViewer);
-PETSC_EXTERN PetscErrorCode DMWritePyLithVertices(DM, PetscViewer);
-PETSC_EXTERN PetscErrorCode DMWritePyLithElements(DM, SectionReal, PetscViewer);
-PETSC_EXTERN PetscErrorCode DMWritePyLithVerticesLocal(DM, PetscViewer);
-PETSC_EXTERN PetscErrorCode DMWritePyLithElementsLocal(DM, SectionReal, PetscViewer);
-
-struct _DMMeshInterpolationInfo {
-  PetscInt   dim;    /*1 The spatial dimension of points */
-  PetscInt   nInput; /* The number of input points */
-  PetscReal *points; /* The input point coordinates */
-  PetscInt  *cells;  /* The cell containing each point */
-  PetscInt   n;      /* The number of local points */
-  Vec        coords; /* The point coordinates */
-  PetscInt   dof;    /* The number of components to interpolate */
-};
-typedef struct _DMMeshInterpolationInfo *DMMeshInterpolationInfo;
-
-PetscErrorCode DMMeshInterpolationCreate(DM dm, DMMeshInterpolationInfo *ctx);
-PetscErrorCode DMMeshInterpolationSetDim(DM dm, PetscInt dim, DMMeshInterpolationInfo ctx);
-PetscErrorCode DMMeshInterpolationGetDim(DM dm, PetscInt *dim, DMMeshInterpolationInfo ctx);
-PetscErrorCode DMMeshInterpolationSetDof(DM dm, PetscInt dof, DMMeshInterpolationInfo ctx);
-PetscErrorCode DMMeshInterpolationGetDof(DM dm, PetscInt *dof, DMMeshInterpolationInfo ctx);
-PetscErrorCode DMMeshInterpolationAddPoints(DM dm, PetscInt n, PetscReal points[], DMMeshInterpolationInfo ctx);
-PetscErrorCode DMMeshInterpolationSetUp(DM dm, DMMeshInterpolationInfo ctx, PetscBool);
-PetscErrorCode DMMeshInterpolationGetCoordinates(DM dm, Vec *points, DMMeshInterpolationInfo ctx);
-PetscErrorCode DMMeshInterpolationGetVector(DM dm, Vec *values, DMMeshInterpolationInfo ctx);
-PetscErrorCode DMMeshInterpolationRestoreVector(DM dm, Vec *values, DMMeshInterpolationInfo ctx);
-PetscErrorCode DMMeshInterpolationEvaluate(DM dm, SectionReal x, Vec v, DMMeshInterpolationInfo ctx);
-PetscErrorCode DMMeshInterpolationDestroy(DM dm, DMMeshInterpolationInfo *ctx);
-
-#endif /* Mesh section */
-#endif

include/petscdmmesh.hh

-#if !defined(__PETSCDMMESH_HH)
-#define __PETSCDMMESH_HH
-
-#include <petscdmmesh.h>
-#include <functional>
-
-using ALE::Obj;
-
-PetscErrorCode DMMeshView_Sieve(const ALE::Obj<PETSC_MESH_TYPE>& mesh, PetscViewer viewer);
-
-#undef __FUNCT__
-#define __FUNCT__ "DMMeshCreateMatrix"
-template<typename Mesh, typename Section>
-PetscErrorCode  DMMeshCreateMatrix(const Obj<Mesh>& mesh, const Obj<Section>& section, MatType mtype, Mat *J, int bs = -1, bool fillMatrix = false)
-{
-  const ALE::Obj<typename Mesh::order_type>& order = mesh->getFactory()->getGlobalOrder(mesh, section->getName(), section);
-  int            localSize  = order->getLocalSize();
-  int            globalSize = order->getGlobalSize();
-  PetscBool      isShell, isBlock, isSeqBlock, isMPIBlock, isSymBlock, isSymSeqBlock, isSymMPIBlock, isSymmetric;
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = MatCreate(mesh->comm(), J);CHKERRQ(ierr);
-  ierr = MatSetSizes(*J, localSize, localSize, globalSize, globalSize);CHKERRQ(ierr);
-  ierr = MatSetType(*J, mtype);CHKERRQ(ierr);
-  ierr = MatSetFromOptions(*J);CHKERRQ(ierr);
-  ierr = PetscStrcmp(mtype, MATSHELL, &isShell);CHKERRQ(ierr);
-  ierr = PetscStrcmp(mtype, MATBAIJ, &isBlock);CHKERRQ(ierr);
-  ierr = PetscStrcmp(mtype, MATSEQBAIJ, &isSeqBlock);CHKERRQ(ierr);
-  ierr = PetscStrcmp(mtype, MATMPIBAIJ, &isMPIBlock);CHKERRQ(ierr);
-  ierr = PetscStrcmp(mtype, MATSBAIJ, &isSymBlock);CHKERRQ(ierr);
-  ierr = PetscStrcmp(mtype, MATSEQSBAIJ, &isSymSeqBlock);CHKERRQ(ierr);
-  ierr = PetscStrcmp(mtype, MATMPISBAIJ, &isSymMPIBlock);CHKERRQ(ierr);
-  // Check for symmetric storage
-  isSymmetric = (PetscBool) (isSymBlock || isSymSeqBlock || isSymMPIBlock);
-  if (!isShell) {
-    PetscInt *dnz, *onz, bsLocal;
-
-    if (bs < 0) {
-      if (isBlock || isSeqBlock || isMPIBlock || isSymBlock || isSymSeqBlock || isSymMPIBlock) {
-        const typename Section::chart_type& chart = section->getChart();
-
-        for(typename Section::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
-          if (section->getFiberDimension(*c_iter)) {
-            bs = section->getFiberDimension(*c_iter);
-            break;
-          }
-        }
-      } else {
-        bs = 1;
-      }
-      // Must have same blocksize on all procs (some might have no points)
-      bsLocal = bs;
-      ierr = MPI_Allreduce(&bsLocal, &bs, 1, MPIU_INT, MPI_MAX, mesh->comm());CHKERRQ(ierr);
-    }
-    ierr = PetscMalloc2(localSize/bs, PetscInt, &dnz, localSize/bs, PetscInt, &onz);CHKERRQ(ierr);
-#ifdef USE_NEW_OVERLAP
-    ierr = preallocateOperatorNewOverlap(mesh, bs, section->getAtlas(), order, dnz, onz, isSymmetric, *J, fillMatrix);CHKERRQ(ierr);
-#else
-    ierr = preallocateOperatorNew(mesh, bs, section->getAtlas(), order, dnz, onz, isSymmetric, *J, fillMatrix);CHKERRQ(ierr);
-#endif
-    ierr = PetscFree2(dnz, onz);CHKERRQ(ierr);
-    if (isSymmetric) {
-      ierr = MatSetOption(*J, MAT_IGNORE_LOWER_TRIANGULAR, PETSC_TRUE);CHKERRQ(ierr);
-    }
-  }
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMMeshCreateGlobalScatter"
-template<typename Mesh, typename Section>
-PetscErrorCode  DMMeshCreateGlobalScatter(const ALE::Obj<Mesh>& m, const ALE::Obj<Section>& s, VecScatter *scatter)
-{
-  const ALE::Obj<typename Mesh::order_type>& globalOrder = m->getFactory()->getGlobalOrder(m, s->getName(), s);
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = DMMeshCreateGlobalScatter(m, s, globalOrder, false, scatter);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMMeshCreateGlobalScatter"
-template<typename Mesh, typename Section>
-PetscErrorCode  DMMeshCreateGlobalScatter(const ALE::Obj<Mesh>& m, const ALE::Obj<Section>& s, const ALE::Obj<typename Mesh::label_type>& label, VecScatter *scatter)
-{
-  const ALE::Obj<typename Mesh::order_type>& globalOrder = m->getFactory()->getGlobalOrder(m, s->getName(), s, -1, label);
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = DMMeshCreateGlobalScatter(m, s, globalOrder, false, scatter);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMMeshCreateGlobalScatter"
-template<typename Mesh, typename Section>
-PetscErrorCode  DMMeshCreateGlobalScatter(const ALE::Obj<Mesh>& m, const std::string& name, const typename Section::chart_type& points, const ALE::Obj<Section>& s, VecScatter *scatter)
-{
-  const ALE::Obj<typename Mesh::order_type>& globalOrder = m->getFactory()->getGlobalOrder(m, name, points, s);
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = DMMeshCreateGlobalScatter(m, s, globalOrder, false, scatter);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "DMMeshCreateGlobalScatter"
-template<typename Mesh, typename Section>
-PetscErrorCode  DMMeshCreateGlobalScatter(const ALE::Obj<Mesh>& m, const ALE::Obj<Section>& s, const ALE::Obj<typename Mesh::order_type>& globalOrder, bool includeConstraints, VecScatter *scatter)
-{
-  typedef typename Mesh::real_section_type::index_type index_type;
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  ierr = PetscLogEventBegin(DMMesh_GetGlobalScatter,0,0,0,0);CHKERRQ(ierr);
-  const typename Mesh::order_type::chart_type& chart = globalOrder->getChart();
-  int *localIndices, *globalIndices;
-  int  localSize   = globalOrder->getLocalSize();
-  int  overlapSize = -1;
-  int  localIndx   = 0, globalIndx = 0;
-  Vec  globalVec, localVec;
-  IS   localIS, globalIS;
-
-  ierr = VecCreate(m->comm(), &globalVec);CHKERRQ(ierr);
-  ierr = VecSetSizes(globalVec, localSize, PETSC_DETERMINE);CHKERRQ(ierr);
-  ierr = VecSetFromOptions(globalVec);CHKERRQ(ierr);
-
-  if (includeConstraints) {
-    overlapSize = s->sizeWithBC();
-    ierr = PetscMalloc(overlapSize*sizeof(int), &localIndices);CHKERRQ(ierr);
-    ierr = PetscMalloc(overlapSize*sizeof(int), &globalIndices);CHKERRQ(ierr);
-  } else {
-    overlapSize = s->size();
-    ierr = PetscMalloc(overlapSize*sizeof(int), &localIndices);CHKERRQ(ierr);
-    ierr = PetscMalloc(overlapSize*sizeof(int), &globalIndices);CHKERRQ(ierr);
-  } // if/else
-
-  // Loop over all local points
-  for(typename Mesh::order_type::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
-    // Map local indices to global indices
-    if (includeConstraints) {
-      s->getIndicesRaw(*p_iter, localIndices, &localIndx, 0);
-      s->getIndicesRaw(*p_iter, globalOrder, globalIndices, &globalIndx, 0);
-    } else {
-      s->getIndices(*p_iter, localIndices, &localIndx, 0, true, true);
-      s->getIndices(*p_iter, globalOrder, globalIndices, &globalIndx, 0, true, false);
-    }
-    //numConstraints += s->getConstraintDimension(*p_iter);
-  }
-  // Local arrays also have constraints, which are not mapped
-  if (localIndx  > overlapSize) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Invalid number of local indices %d, should not be greater than %d", localIndx, overlapSize);
-  if (globalIndx > overlapSize) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Invalid number of global indices %d, should not be greater than %d", globalIndx, overlapSize);
-  if (globalIndx != localIndx)  SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Mismatched number of global indices %d, and local indices %d", globalIndx, localIndx);
-  if (m->debug()) {
-    globalOrder->view("Global Order");
-    for(int i = 0; i < globalIndx; ++i) {
-      printf("[%d] localIndex[%d]: %d globalIndex[%d]: %d\n", m->commRank(), i, localIndices[i], i, globalIndices[i]);
-    }
-  }
-  ierr = ISCreateGeneral(PETSC_COMM_SELF, localIndx, localIndices,PETSC_OWN_POINTER,  &localIS);CHKERRQ(ierr);
-  ierr = ISCreateGeneral(PETSC_COMM_SELF, globalIndx, globalIndices,PETSC_OWN_POINTER, &globalIS);CHKERRQ(ierr);
-  // Can remove this when I test it with NULL
-#ifdef PETSC_USE_COMPLEX
-  ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, 1,s->getStorageSize(), NULL, &localVec);CHKERRQ(ierr);
-#else
-  ierr = VecCreateSeqWithArray(PETSC_COMM_SELF, 1,s->getStorageSize(), s->restrictSpace(), &localVec);CHKERRQ(ierr);
-#endif
-  ierr = VecScatterCreate(localVec, localIS, globalVec, globalIS, scatter);CHKERRQ(ierr);
-  ierr = ISDestroy(&globalIS);CHKERRQ(ierr);
-  ierr = ISDestroy(&localIS);CHKERRQ(ierr);
-  ierr = VecDestroy(&localVec);CHKERRQ(ierr);
-  ierr = VecDestroy(&globalVec);CHKERRQ(ierr);
-  ierr = PetscLogEventEnd(DMMesh_GetGlobalScatter,0,0,0,0);CHKERRQ(ierr);
-  PetscFunctionReturn(0);
-}
-
-template<typename Mesh, typename Section>
-void createOperator(const ALE::Obj<Mesh>& mesh, const ALE::Obj<Section>& s, const ALE::Obj<Mesh>& op) {
-  typedef ALE::SieveAlg<Mesh> sieve_alg_type;
-  typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
-  const typename Section::chart_type& chart = s->getChart();
-
-  // Create local operator
-  //   We do not decorate arrows yet
-  for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
-    const Obj<typename sieve_alg_type::supportArray>& star = sieve_alg_type::star(mesh, *p_iter);
-
-    for(typename sieve_alg_type::supportArray::const_iterator s_iter = star->begin(); s_iter != star->end(); ++s_iter) {
-      const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(mesh, *s_iter);
-
-      for(typename sieve_alg_type::coneArray::const_iterator c_iter = closure->begin(); c_iter != closure->end(); ++c_iter) {
-        op->getSieve()->addCone(*c_iter, *p_iter);
-      }
-    }
-  }
-  op->view("Local operator");
-  // Construct overlap
-  Obj<FlexMesh::send_overlap_type> sendOverlap = mesh->getSendOverlap();
-  Obj<FlexMesh::recv_overlap_type> recvOverlap = mesh->getRecvOverlap();
-  FlexMesh::renumbering_type&      renumbering = mesh->getRenumbering();
-
-  sendOverlap->view("Mesh send overlap");
-  recvOverlap->view("Mesh recv overlap");
-  // Complete operator
-  typedef ALE::DistributionNew<FlexMesh>::cones_type ConeOverlap;
-
-  ALE::Obj<ConeOverlap> overlapCones = ALE::DistributionNew<FlexMesh>::completeCones(op->getSieve(), op->getSieve(), renumbering, sendOverlap, recvOverlap);
-  op->view("Completed operator");
-  // Update renumbering and overlap
-  overlapCones->view("Overlap cones");
-  Obj<FlexMesh::send_overlap_type>        opSendOverlap = op->getSendOverlap();
-  Obj<FlexMesh::recv_overlap_type>        opRecvOverlap = op->getRecvOverlap();
-  FlexMesh::renumbering_type&             opRenumbering = op->getRenumbering();
-  const typename ConeOverlap::chart_type& overlapChart  = overlapCones->getChart();
-  int                                     p             = renumbering.size();
-
-  opRenumbering = renumbering;
-  for(typename ConeOverlap::chart_type::const_iterator p_iter = overlapChart.begin(); p_iter != overlapChart.end(); ++p_iter) {
-    if (opRenumbering.find(p_iter->second) == opRenumbering.end()) {
-      opRenumbering[p_iter->second] = p++;
-    }
-  }
-  ALE::SetFromMap<FlexMesh::renumbering_type> opGlobalPoints(opRenumbering);
-
-  ALE::OverlapBuilder<>::constructOverlap(opGlobalPoints, opRenumbering, opSendOverlap, opRecvOverlap);
-  sendOverlap->view("Operator send overlap");
-  recvOverlap->view("Operator recv overlap");
-  // Create global order
-  Obj<FlexMesh::order_type> globalOrder = new FlexMesh::order_type(op->comm(), op->debug());
-
-  op->getFactory()->constructLocalOrder(globalOrder, opSendOverlap, opGlobalPoints, s);
-  op->getFactory()->calculateOffsets(globalOrder);
-  op->getFactory()->updateOrder(globalOrder, opGlobalPoints);
-  op->getFactory()->completeOrder(globalOrder, opSendOverlap, opRecvOverlap);
-  globalOrder->view("Operator global order");
-  // Create dnz/onz or CSR
-};
-
-template<typename Atlas>
-class AdjVisitor {
-public:
-  typedef typename ALE::Mesh<PetscInt,PetscScalar>::point_type point_type;
-protected:
-  Atlas&                 atlas;
-  ALE::Mesh<PetscInt,PetscScalar>::sieve_type& adjGraph;
-  point_type             p;
-public:
-  AdjVisitor(Atlas& atlas, ALE::Mesh<PetscInt,PetscScalar>::sieve_type& adjGraph) : atlas(atlas), adjGraph(adjGraph) {};
-  void visitPoint(const point_type& point) {
-    if (atlas.restrictPoint(point)[0].prefix) {
-      adjGraph.addCone(point, p);
-    }
-  };
-  template<typename Arrow>
-  void visitArrow(const Arrow&) {};
-public:
-  void setRoot(const point_type& point) {this->p = point;};
-};
-
-#undef __FUNCT__
-#define __FUNCT__ "createLocalAdjacencyGraph"
-template<typename Mesh, typename Atlas>
-PetscErrorCode createLocalAdjacencyGraph(const ALE::Obj<Mesh>& mesh, const ALE::Obj<Atlas>& atlas, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph)
-{
-  typedef typename ALE::ISieveVisitor::TransitiveClosureVisitor<typename Mesh::sieve_type, AdjVisitor<Atlas> > ClosureVisitor;
-  typedef typename ALE::ISieveVisitor::ConeVisitor<typename Mesh::sieve_type, ClosureVisitor>                  ConeVisitor;
-  typedef typename ALE::ISieveVisitor::TransitiveClosureVisitor<typename Mesh::sieve_type, ConeVisitor>        StarVisitor;
-  AdjVisitor<Atlas> adjV(*atlas, *adjGraph);
-  ClosureVisitor    closureV(*mesh->getSieve(), adjV);
-  ConeVisitor       coneV(*mesh->getSieve(), closureV);
-  StarVisitor       starV(*mesh->getSieve(), coneV);
-  /* In general, we need to get FIAT info that attaches dual basis vectors to sieve points */
-  const typename Atlas::chart_type& chart = atlas->getChart();
-
-  PetscFunctionBegin;
-  starV.setIsCone(false);
-  for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
-    adjV.setRoot(*c_iter);
-    mesh->getSieve()->support(*c_iter, starV);
-    closureV.clear();
-    starV.clear();
-  }
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "createLocalAdjacencyGraphI"
-template<typename Mesh, typename Atlas>
-PetscErrorCode createLocalAdjacencyGraphI(const ALE::Obj<Mesh>& mesh, const ALE::Obj<Atlas>& atlas, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph)
-{
-  typedef typename ALE::ISieveVisitor::TransitiveClosureVisitor<typename Mesh::sieve_type, AdjVisitor<Atlas> > ClosureVisitor;
-  typedef typename ALE::ISieveVisitor::ConeVisitor<typename Mesh::sieve_type, ClosureVisitor>                  ConeVisitor;
-  typedef typename ALE::ISieveVisitor::TransitiveClosureVisitor<typename Mesh::sieve_type, ConeVisitor>        StarVisitor;
-  AdjVisitor<Atlas> adjV(*atlas, *adjGraph);
-  ClosureVisitor    closureV(*mesh->getSieve(), adjV);
-  ConeVisitor       coneV(*mesh->getSieve(), closureV);
-  StarVisitor       starV(*mesh->getSieve(), coneV);
-  /* In general, we need to get FIAT info that attaches dual basis vectors to sieve points */
-  const typename Atlas::chart_type& chart = atlas->getChart();
-
-  PetscFunctionBegin;
-  // Changes for ISieve
-  //   1) Add AdjSizeVisitor to set cone sizes
-  //   2) Add new symmetrizeSizes() to ISieve to calculate support sizes
-  //   3) Allocate adjGraph
-  //   4) Change AdjVisitor to stack up cone rather than calling addPoint()
-  //   5) Get points and call setCone() each time
-  starV.setIsCone(false);
-  for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
-    adjV.setRoot(*c_iter);
-    mesh->getSieve()->support(*c_iter, starV);
-    closureV.clear();
-    starV.clear();
-  }
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "createLocalAdjacencyGraph"
-template<typename Atlas>
-PetscErrorCode createLocalAdjacencyGraph(const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar> >& mesh, const ALE::Obj<Atlas>& atlas, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph)
-{
-  typedef ALE::SieveAlg<ALE::Mesh<PetscInt,PetscScalar> > sieve_alg_type;
-  /* In general, we need to get FIAT info that attaches dual basis vectors to sieve points */
-  const typename Atlas::chart_type& chart = atlas->getChart();
-
-  PetscFunctionBegin;
-  for(typename Atlas::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
-    const Obj<typename sieve_alg_type::supportArray>& star = sieve_alg_type::star(mesh, *c_iter);
-
-    for(typename sieve_alg_type::supportArray::const_iterator s_iter = star->begin(); s_iter != star->end(); ++s_iter) {
-      const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(mesh, *s_iter);
-
-      for(typename sieve_alg_type::coneArray::const_iterator cl_iter = closure->begin(); cl_iter != closure->end(); ++cl_iter) {
-        adjGraph->addCone(*cl_iter, *c_iter);
-      }
-    }
-  }
-  PetscFunctionReturn(0);
-}
-
-template<typename Arrow>
-struct SelectSource : public std::unary_function<Arrow, typename Arrow::source_type>
-{
-  typename Arrow::source_type& operator()(Arrow& x) const {return x.source;}
-  const typename Arrow::source_type& operator()(const Arrow& x) const {return x.source;}
-};
-
-template<class Pair>
-struct Select1st : public std::unary_function<Pair, typename Pair::first_type>
-{
-  typename Pair::first_type& operator()(Pair& x) const {return x.first;}
-  const typename Pair::first_type& operator()(const Pair& x) const {return x.first;}
-};
-
-template<typename Mesh, typename Order>
-PetscErrorCode globalizeLocalAdjacencyGraph(const ALE::Obj<Mesh>& mesh, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::send_overlap_type>& sendOverlap, const ALE::Obj<Order>& globalOrder)
-{
-  typedef ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
-  const int debug = mesh->debug();
-
-  PetscFunctionBegin;
-  ALE::PointFactory<typename Mesh::point_type>& pointFactory = ALE::PointFactory<FlexMesh::point_type>::singleton(mesh->comm(), mesh->getSieve()->getChart().max(), mesh->debug());
-  // Check for points not in sendOverlap
-  std::set<typename Mesh::point_type> interiorPoints;
-  std::set<typename Mesh::point_type> overlapSources;
-  std::set<typename Mesh::sieve_type::arrow_type> overlapArrows;
-  const Obj<FlexMesh::sieve_type::traits::capSequence>& columns = adjGraph->cap();
-
-  for(FlexMesh::sieve_type::traits::capSequence::iterator p_iter = columns->begin(); p_iter != columns->end(); ++p_iter) {
-    // This optimization does not work because empty points are sent anyway
-    //if (!sendOverlap->support(*p_iter)->size() && globalOrder->restrictPoint(*p_iter)[0].index) {
-    if (!sendOverlap->support(*p_iter)->size()) {
-      interiorPoints.insert(*p_iter);
-    } else {
-      // If a point p is in the overlap for process i and an adjacent point q is in the overlap for process j then:
-      //   Replace (q, p) with (q', p)
-      //   Notice I can reverse the arrow because the graph is initially symmetric
-      if (debug) {std::cout << "["<<globalOrder->commRank()<<"] Checking overlap point " << *p_iter << " for overlap neighbors" << std::endl;}
-      const Obj<typename FlexMesh::sieve_type::supportSequence>&     neighbors = adjGraph->support(*p_iter);
-      const typename FlexMesh::sieve_type::supportSequence::iterator nEnd      = neighbors->end();
-
-      for(typename FlexMesh::sieve_type::supportSequence::iterator n_iter = neighbors->begin(); n_iter != nEnd; ++n_iter) {
-        if (sendOverlap->support(*n_iter)->size()) {
-          if (debug) {std::cout << "["<<globalOrder->commRank()<<"]   Found overlap neighbor " << *n_iter << std::endl;}
-          const Obj<typename FlexMesh::send_overlap_type::supportSequence>&     ranks = sendOverlap->support(*p_iter);
-          const typename FlexMesh::send_overlap_type::supportSequence::iterator rEnd  = ranks->end();
-          bool equal = true;
-
-          for(typename FlexMesh::send_overlap_type::supportSequence::iterator r_iter = ranks->begin(); r_iter != rEnd; ++r_iter) {
-            const Obj<typename FlexMesh::send_overlap_type::supportSequence>&     nRanks = sendOverlap->support(*n_iter);
-            const typename FlexMesh::send_overlap_type::supportSequence::iterator nrEnd  = nRanks->end();
-            bool found = false;
-
-            if (debug) {std::cout << "["<<globalOrder->commRank()<<"]     Checking overlap rank " << *r_iter << std::endl;}
-            for(typename FlexMesh::send_overlap_type::supportSequence::iterator nr_iter = nRanks->begin(); nr_iter != nrEnd; ++nr_iter) {
-              if (debug) {std::cout << "["<<globalOrder->commRank()<<"]     Checking neighbor overlap rank " << *nr_iter << std::endl;}
-              if (*nr_iter == *r_iter) {
-                found = true;
-                break;
-              }
-            }
-            if (!found) {
-              equal = false;
-              break;
-            }
-          }
-          if (!equal) {
-            if (debug) {std::cout << "["<<globalOrder->commRank()<<"]   Unequal rank sets, replacing arrow " << *n_iter <<" --> "<<*p_iter << std::endl;}
-            overlapArrows.insert(typename Mesh::sieve_type::arrow_type(*n_iter, *p_iter));
-          } else {
-            if (debug) {std::cout << "["<<globalOrder->commRank()<<"]   Equal rank sets, doing nothing for arrow " << *n_iter <<" --> "<<*p_iter << std::endl;}
-          }
-        }
-      }
-    }
-  }
-  // Renumber those points
-  pointFactory.clear();
-  pointFactory.setMax(mesh->getSieve()->getChart().max());
-  pointFactory.renumberPoints(interiorPoints.begin(), interiorPoints.end());
-  //pointFactory.renumberPoints(overlapArrows.begin(), overlapArrows.end(), SelectSource<typename Mesh::sieve_type::arrow_type>());
-  // They should use a key extractor: overlapSources.insert(overlapArrows.begin(), overlapArrows.end(), SelectSource<typename Mesh::sieve_type::arrow_type>());
-  for(typename std::set<typename Mesh::sieve_type::arrow_type>::const_iterator a_iter = overlapArrows.begin(); a_iter != overlapArrows.end(); ++a_iter) {
-    overlapSources.insert(a_iter->source);
-  }
-  pointFactory.renumberPoints(overlapSources.begin(), overlapSources.end());
-  typename ALE::PointFactory<typename Mesh::point_type>::renumbering_type& renumbering    = pointFactory.getRenumbering();
-  typename ALE::PointFactory<typename Mesh::point_type>::renumbering_type& invRenumbering = pointFactory.getInvRenumbering();
-  // Replace points in local adjacency graph
-  const Obj<FlexMesh::sieve_type::traits::baseSequence>& base    = adjGraph->base();
-  ALE::Obj<std::vector<typename Mesh::point_type> >       newCone = new std::vector<typename Mesh::point_type>();
-
-  for(FlexMesh::sieve_type::traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
-    const ALE::Obj<FlexMesh::sieve_type::coneSequence>& cone    = adjGraph->cone(*b_iter);
-    const FlexMesh::sieve_type::coneSequence::iterator  cEnd    = cone->end();
-    bool                                                 replace = false;
-
-    for(FlexMesh::sieve_type::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
-      if (interiorPoints.find(*c_iter) != interiorPoints.end()) {
-        newCone->push_back(invRenumbering[*c_iter]);
-        replace = true;
-      } else {
-        newCone->push_back(*c_iter);
-      }
-    }
-    if (interiorPoints.find(*b_iter) != interiorPoints.end()) {
-      adjGraph->clearCone(*b_iter);
-      adjGraph->setCone(newCone, invRenumbering[*b_iter]);
-      if (debug) {std::cout << "["<<globalOrder->commRank()<<"]: Replacing cone for " << *b_iter << "("<<invRenumbering[*b_iter]<<") with" << std::endl;}
-    } else if (replace) {
-      adjGraph->clearCone(*b_iter);
-      adjGraph->setCone(newCone, *b_iter);
-      if (debug) {std::cout << "["<<globalOrder->commRank()<<"]: Replacing cone for " << *b_iter << " with" << std::endl;}
-    }
-    if (debug && ((interiorPoints.find(*b_iter) != interiorPoints.end()) || replace)) {
-      for(typename std::vector<typename Mesh::point_type>::const_iterator c_iter = newCone->begin(); c_iter != newCone->end(); ++c_iter) {
-        std::cout << "["<<globalOrder->commRank()<<"]:   point " << *c_iter << std::endl;
-      }
-    }
-    newCone->clear();
-  }
-  // Replace arrows
-  for(typename std::set<typename Mesh::sieve_type::arrow_type>::const_iterator a_iter = overlapArrows.begin(); a_iter != overlapArrows.end(); ++a_iter) {
-    if (debug) {std::cout << "["<<globalOrder->commRank()<<"]: Replacing " << a_iter->source<<" --> "<<a_iter->target << " with " << invRenumbering[a_iter->source]<<" --> "<<a_iter->target << std::endl;}
-    adjGraph->removeArrow(a_iter->source, a_iter->target);
-    adjGraph->addArrow(invRenumbering[a_iter->source], a_iter->target);
-  }
-  // Add new points into ordering
-#if 1
-  for(typename ALE::PointFactory<typename Mesh::point_type>::renumbering_type::const_iterator p_iter = renumbering.begin(); p_iter != renumbering.end(); ++p_iter) {
-    if (debug) {std::cout << "["<<globalOrder->commRank()<<"]: Updating " << p_iter->first << " to " << globalOrder->restrictPoint(p_iter->second)[0] << std::endl;}
-    globalOrder->addPoint(p_iter->first);
-    globalOrder->updatePoint(p_iter->first, globalOrder->restrictPoint(p_iter->second));
-  }
-#else
-  globalOrder->reallocatePoint(renumbering.begin(), renumbering.end(), Select1st<typename ALE::PointFactory<typename Mesh::point_type>::renumbering_type::const_iterator::value_type>());
-  for(typename ALE::PointFactory<typename Mesh::point_type>::renumbering_type::const_iterator p_iter = renumbering.begin(); p_iter != renumbering.end(); ++p_iter) {
-    if (debug) {std::cout << "["<<globalOrder->commRank()<<"]: Updating " << p_iter->first << " to " << globalOrder->restrictPoint(p_iter->second)[0] << std::endl;}
-    globalOrder->updatePoint(p_iter->first, globalOrder->restrictPoint(p_iter->second));
-  }
-#endif
-  PetscFunctionReturn(0);
-}
-
-template<typename Order>
-PetscErrorCode globalizeLocalAdjacencyGraph(const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar> >& mesh, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::send_overlap_type>& sendOverlap, const ALE::Obj<Order>& globalOrder)
-{
-  PetscFunctionBegin;
-  PetscFunctionReturn(0);
-}
-
-template<typename Mesh, typename Order>
-PetscErrorCode localizeLocalAdjacencyGraph(const ALE::Obj<Mesh>& mesh, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::send_overlap_type>& sendOverlap, const ALE::Obj<Order>& globalOrder)
-{
-  PetscFunctionBegin;
-  ALE::PointFactory<typename Mesh::point_type>& pointFactory = ALE::PointFactory<ALE::Mesh<PetscInt,PetscScalar>::point_type>::singleton(mesh->comm(), mesh->getSieve()->getChart().max(), mesh->debug());
-  typename ALE::PointFactory<typename Mesh::point_type>::renumbering_type& renumbering = pointFactory.getRenumbering();
-  // Replace points in local adjacency graph
-  const Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type::traits::baseSequence>& base = adjGraph->base();
-
-  for(ALE::Mesh<PetscInt,PetscScalar>::sieve_type::traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
-    // Replace globalized points
-    if (renumbering.find(*b_iter) != renumbering.end()) {
-      adjGraph->clearCone(renumbering[*b_iter]);
-      adjGraph->setCone(adjGraph->cone(*b_iter), renumbering[*b_iter]);
-      adjGraph->clearCone(*b_iter);
-    }
-  }
-  PetscFunctionReturn(0);
-}
-
-template<typename Order>
-PetscErrorCode localizeLocalAdjacencyGraph(const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar> >& mesh, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::send_overlap_type>& sendOverlap, const ALE::Obj<Order>& globalOrder)
-{
-  PetscFunctionBegin;
-  PetscFunctionReturn(0);
-}
-
-template<typename Mesh, typename Order>
-PetscErrorCode renumberLocalAdjacencyGraph(const ALE::Obj<Mesh>& mesh, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::send_overlap_type>& sendOverlap, const ALE::Obj<Order>& globalOrder)
-{
-  typedef typename ALE::Mesh<PetscInt,PetscScalar> FlexMesh;
-  ALE::Obj<std::set<typename Mesh::point_type> > newCone = new std::set<typename Mesh::point_type>();
-  const int debug = mesh->debug();
-
-  PetscFunctionBegin;
-  ALE::PointFactory<typename Mesh::point_type>& pointFactory = ALE::PointFactory<typename FlexMesh::point_type>::singleton(mesh->comm(), mesh->getSieve()->getChart().max(), debug);
-  pointFactory.constructRemoteRenumbering();
-  typename ALE::PointFactory<typename Mesh::point_type>::renumbering_type&        renumbering       = pointFactory.getRenumbering();
-  typename ALE::PointFactory<typename Mesh::point_type>::remote_renumbering_type& remoteRenumbering = pointFactory.getRemoteRenumbering();
-  // Replace points in local adjacency graph
-  const Obj<typename FlexMesh::sieve_type::traits::baseSequence>& base = adjGraph->base();
-
-  for(FlexMesh::sieve_type::traits::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
-    // Loop over cone checking for remote points that shadow local points
-    const Obj<FlexMesh::sieve_type::traits::coneSequence>&     cone = adjGraph->cone(*b_iter);
-    const FlexMesh::sieve_type::traits::coneSequence::iterator cEnd = cone->end();
-    bool replace = false;
-
-    if (debug) {std::cout <<"["<<adjGraph->commRank()<<"]: Checking base point " << *b_iter << std::endl;}
-    for(FlexMesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
-      bool useOldPoint = true;
-
-      if (debug) {std::cout <<"["<<adjGraph->commRank()<<"]:   cone point " << *c_iter;}
-      if (!globalOrder->isLocal(*c_iter)) {
-        if (debug) {std::cout << " is nonlocal" << std::endl;}
-        for(int p = 0; p < adjGraph->commSize(); ++p) {
-          if (p == adjGraph->commRank()) continue;
-          if (remoteRenumbering[p].find(*c_iter) != remoteRenumbering[p].end()) {
-            if (debug) {std::cout <<"["<<adjGraph->commRank()<<"]:   found " << *c_iter << " on process " << p << " as point " << remoteRenumbering[p][*c_iter];}
-            const Obj<FlexMesh::send_overlap_type::traits::coneSequence>& localPoint = sendOverlap->cone(p, remoteRenumbering[p][*c_iter]);
-
-            if (localPoint->size()) {
-              if (debug) {std::cout << " with local match " << *localPoint->begin() << std::endl;}
-              newCone->insert(*localPoint->begin());
-              replace     = true;
-              useOldPoint = false;
-              break;
-            } else {
-              if (debug) {std::cout << " but not in sendOverlap" << std::endl;}
-            }
-          }
-        }
-      } else {
-        if (debug) {std::cout << " is local" << std::endl;}
-        if (renumbering.find(*c_iter) != renumbering.end()) {
-          if (debug) {std::cout <<"["<<adjGraph->commRank()<<"]:   found " << *c_iter << " locally as point " << renumbering[*c_iter];}
-          newCone->insert(renumbering[*c_iter]);
-          replace     = true;
-          useOldPoint = false;
-        }
-      }
-      if (useOldPoint) {
-        newCone->insert(*c_iter);
-      }
-    }
-    if (replace) {
-      if (debug > 1) {
-        std::cout <<"["<<adjGraph->commRank()<<"]: Replacing cone for " << *b_iter << " due to shadowed points from" << std::endl;
-        const Obj<FlexMesh::sieve_type::traits::coneSequence>&     cone = adjGraph->cone(*b_iter);
-        const FlexMesh::sieve_type::traits::coneSequence::iterator cEnd = cone->end();
-        for(typename FlexMesh::sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cEnd; ++c_iter) {
-          std::cout <<"["<<adjGraph->commRank()<<"]:   point " << *c_iter << std::endl;
-        }
-        std::cout <<"["<<adjGraph->commRank()<<"]: to" << std::endl;
-        for(typename std::set<typename Mesh::point_type>::const_iterator c_iter = newCone->begin(); c_iter != newCone->end(); ++c_iter) {
-          std::cout <<"["<<adjGraph->commRank()<<"]:   point " << *c_iter << std::endl;
-        }
-      }
-      adjGraph->setCone(newCone, *b_iter);
-      newCone->clear();
-    }
-  }
-  PetscFunctionReturn(0);
-}
-
-template<typename Order>
-PetscErrorCode renumberLocalAdjacencyGraph(const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar> >& mesh, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type>& adjGraph, const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::send_overlap_type>& sendOverlap, const ALE::Obj<Order>& globalOrder)
-{
-  PetscFunctionBegin;
-  PetscFunctionReturn(0);
-}
-
-#undef __FUNCT__
-#define __FUNCT__ "preallocateOperator"
-/* Problem:
-     We want to allocate an operator. This means we want to know the ordering of all unknowns in the sparsity pattern.
-     The preexisting overlap will not contain information about all unknowns (columns) in the operator.
-
-   Solution:
-     Construct the local sparsity pattern, using globally consistent names for new interior points. Cone complete this
-     pattern, which gives an augmented overlap structure. Insert offsets for the new, global point ids in the existing
-     order, and then complete the order. This argues for a recreation of the order, rather than use of an existing
-     order.
-
-   Problem:
-     Figure out sparsity pattern of the operator, when we have already locally numbered all points. The overlap can
-     establish common names for shared points, but not for interior points.
-
-   Solution:
-     Create a shared resource that bestows globally consistent names.
-*/
-template<typename Mesh, typename Atlas>
-PetscErrorCode preallocateOperator(const ALE::Obj<Mesh>& mesh, const int bs, const ALE::Obj<Atlas>& atlas, const ALE::Obj<typename Mesh::order_type>& globalOrder, PetscInt dnz[], PetscInt onz[], Mat A)
-{
-  MPI_Comm                              comm      = mesh->comm();
-  const int                             rank      = mesh->commRank();
-  const int                             debug     = mesh->debug();
-  const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar> >            adjBundle = new ALE::Mesh<PetscInt,PetscScalar>(comm, debug);
-  const ALE::Obj<ALE::Mesh<PetscInt,PetscScalar>::sieve_type> adjGraph  = new ALE::Mesh<PetscInt,PetscScalar>::sieve_type(comm, debug);
-  PetscInt       numLocalRows, firstRow;
-  PetscErrorCode ierr;
-
-  PetscFunctionBegin;
-  adjBundle->setSieve(adjGraph);
-  numLocalRows = globalOrder->getLocalSize();
-  firstRow     = globalOrder->getGlobalOffsets()[rank];
-  ierr = createLocalAdjacencyGraph(mesh, atlas, adjGraph);CHKERRQ(ierr);
-  /* Distribute adjacency graph */
-  adjBundle->constructOverlap();
-  typedef typename Mesh::sieve_type::point_type point_type;
-  typedef typename Mesh::send_overlap_type      send_overlap_type;
-  typedef typename Mesh::recv_overlap_type      recv_overlap_type;
-  typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, point_type> > send_section_type;
-  typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, point_type> > recv_section_type;
-  const Obj<send_overlap_type>& vertexSendOverlap = mesh->getSendOverlap();
-  const Obj<recv_overlap_type>& vertexRecvOverlap = mesh->getRecvOverlap();
-  const Obj<send_overlap_type>  nbrSendOverlap    = new send_overlap_type(comm, debug);
-  const Obj<recv_overlap_type>  nbrRecvOverlap    = new recv_overlap_type(comm, debug);
-  const Obj<send_section_type>  sendSection       = new send_section_type(comm, debug);
-  const Obj<recv_section_type>  recvSection       = new recv_section_type(comm, sendSection->getTag(), debug);
-