I'm not sure exactly what you mean, but I'll try to give more details. We have our own DM class (DM_Moose) and we set our own field and domain decomposition routines:
dm->ops->createfielddecomposition = DMCreateFieldDecomposition_Moose;
dm->ops->createdomaindecomposition = DMCreateDomainDecomposition_Moose;
static PetscErrorCode
DMCreateFieldDecomposition_Moose(
DM dm, PetscInt * len, char *** namelist, IS ** islist, DM ** dmlist)
{
PetscErrorCode ierr;
DM_Moose * dmm = (DM_Moose *)(dm->data);
PetscFunctionBegin;
/* Only called after DMSetUp(). */
if (!dmm->_splitlocs)
PetscFunctionReturn(0);
*len = dmm->_splitlocs->size();
if (namelist)
{
ierr = PetscMalloc(*len * sizeof(char *), namelist);
CHKERRQ(ierr);
}
if (islist)
{
ierr = PetscMalloc(*len * sizeof(IS), islist);
CHKERRQ(ierr);
}
if (dmlist)
{
ierr = PetscMalloc(*len * sizeof(DM), dmlist);
CHKERRQ(ierr);
}
for (const auto & dit : *(dmm->_splitlocs))
{
unsigned int d = dit.second;
std::string dname = dit.first;
DM_Moose::SplitInfo & dinfo = (*dmm->_splits)[dname];
if (!dinfo._dm)
{
ierr = DMCreateMoose(((PetscObject)dm)->comm, *dmm->_nl, &dinfo._dm);
CHKERRQ(ierr);
ierr = PetscObjectSetOptionsPrefix((PetscObject)dinfo._dm, ((PetscObject)dm)->prefix);
CHKERRQ(ierr);
std::string suffix = std::string("fieldsplit_") + dname + "_";
ierr = PetscObjectAppendOptionsPrefix((PetscObject)dinfo._dm, suffix.c_str());
CHKERRQ(ierr);
}
ierr = DMSetFromOptions(dinfo._dm);
CHKERRQ(ierr);
ierr = DMSetUp(dinfo._dm);
CHKERRQ(ierr);
if (namelist)
{
ierr = PetscStrallocpy(dname.c_str(), (*namelist) + d);
CHKERRQ(ierr);
}
if (islist)
{
if (!dinfo._rembedding)
{
IS dembedding, lembedding;
ierr = DMMooseGetEmbedding_Private(dinfo._dm, &dembedding);
CHKERRQ(ierr);
if (dmm->_embedding)
{
// Create a relative embedding into the parent's index space.
ierr = ISEmbed(dembedding, dmm->_embedding, PETSC_TRUE, &lembedding);
CHKERRQ(ierr);
const PetscInt * lindices;
PetscInt len, dlen, llen, *rindices, off, i;
ierr = ISGetLocalSize(dembedding, &dlen);
CHKERRQ(ierr);
ierr = ISGetLocalSize(lembedding, &llen);
CHKERRQ(ierr);
if (llen != dlen)
SETERRQ1(((PetscObject)dm)->comm, PETSC_ERR_PLIB, "Failed to embed split %D", d);
ierr = ISDestroy(&dembedding);
CHKERRQ(ierr);
// Convert local embedding to global (but still relative) embedding
ierr = PetscMalloc(llen * sizeof(PetscInt), &rindices);
CHKERRQ(ierr);
ierr = ISGetIndices(lembedding, &lindices);
CHKERRQ(ierr);
ierr = PetscMemcpy(rindices, lindices, llen * sizeof(PetscInt));
CHKERRQ(ierr);
ierr = ISDestroy(&lembedding);
CHKERRQ(ierr);
// We could get the index offset from a corresponding global vector, but subDMs don't yet
// have global vectors
ierr = ISGetLocalSize(dmm->_embedding, &len);
CHKERRQ(ierr);
ierr = MPI_Scan(&len,
&off,
1,
#ifdef PETSC_USE_64BIT_INDICES
MPI_LONG_LONG_INT,
#else
MPI_INT,
#endif
MPI_SUM,
((PetscObject)dm)->comm);
CHKERRQ(ierr);
off -= len;
for (i = 0; i < llen; ++i)
rindices[i] += off;
ierr = ISCreateGeneral(
((PetscObject)dm)->comm, llen, rindices, PETSC_OWN_POINTER, &(dinfo._rembedding));
CHKERRQ(ierr);
}
else
{
dinfo._rembedding = dembedding;
}
}
ierr = PetscObjectReference((PetscObject)(dinfo._rembedding));
CHKERRQ(ierr);
(*islist)[d] = dinfo._rembedding;
}
if (dmlist)
{
ierr = PetscObjectReference((PetscObject)dinfo._dm);
CHKERRQ(ierr);
(*dmlist)[d] = dinfo._dm;
}
}
PetscFunctionReturn(0);
}
static PetscErrorCode
DMCreateDomainDecomposition_Moose(
DM dm, PetscInt * len, char *** namelist, IS ** innerislist, IS ** outerislist, DM ** dmlist)
{
PetscErrorCode ierr;
PetscFunctionBegin;
/* Use DMCreateFieldDecomposition_Moose() to obtain everything but outerislist, which is currently
* PETSC_NULL. */
if (outerislist)
*outerislist = PETSC_NULL; /* FIX: allow mesh-based overlap. */
ierr = DMCreateFieldDecomposition_Moose(dm, len, namelist, innerislist, dmlist);
CHKERRQ(ierr);
PetscFunctionReturn(0);
}