Actual source code: dmmoab.cxx

      1: #include <petsc/private/dmmbimpl.h>

      3: #include <petscdmmoab.h>
      4: #include <MBTagConventions.hpp>
      5: #include <moab/NestedRefine.hpp>
      6: #include <moab/Skinner.hpp>

      8: /*MC
      9:   DMMOAB = "moab" - A `DM` object that encapsulates an unstructured mesh described by the MOAB mesh database {cite}`moabwebsite`.
     10:                     Direct access to the MOAB Interface and other mesh manipulation related objects are available
     11:                     through public API. Ability to create global and local representation of Vecs containing all
     12:                     unknowns in the interior and shared boundary via a transparent tag-data wrapper is provided
     13:                     along with utility functions to traverse the mesh and assemble a discrete system via
     14:                     field-based/blocked Vec(Get/Set) methods. Input from and output to different formats are
     15:                     available.

     17:   Level: intermediate

     19: .seealso: `DMMOAB`, `DMType`, `DMMoabCreate()`, `DMCreate()`, `DMSetType()`, `DMMoabCreateMoab()`
     20: M*/

     22: /* External function declarations here */
     23: PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
     24: PETSC_EXTERN PetscErrorCode DMCreateDefaultConstraints_Moab(DM dm);
     25: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
     26: PETSC_EXTERN PetscErrorCode DMCreateCoordinateDM_Moab(DM dm, DM *cdm);
     27: PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM *dmRefined);
     28: PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM *dmCoarsened);
     29: PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmRefined[]);
     30: PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmCoarsened[]);
     31: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm);
     32: PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM, Vec *);
     33: PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM, Vec *);
     34: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
     35: PETSC_EXTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM, Vec, InsertMode, Vec);
     36: PETSC_EXTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM, Vec, InsertMode, Vec);
     37: PETSC_EXTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM, Vec, InsertMode, Vec);
     38: PETSC_EXTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM, Vec, InsertMode, Vec);

     40: /* Un-implemented routines */
     41: /*
     42: PETSC_EXTERN PetscErrorCode DMCreatelocalsection_Moab(DM dm);
     43: PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dmCoarse, DM dmFine, Mat *mat);
     44: PETSC_EXTERN PetscErrorCode DMLoad_Moab(DM dm, PetscViewer viewer);
     45: PETSC_EXTERN PetscErrorCode DMGetDimPoints_Moab(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd);
     46: PETSC_EXTERN PetscErrorCode DMCreateSubDM_Moab(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
     47: PETSC_EXTERN PetscErrorCode DMLocatePoints_Moab(DM dm, Vec v, IS *cellIS);
     48: */

     50: /*@C
     51:   DMMoabCreate - Creates a `DMMOAB` object, which encapsulates a moab instance

     53:   Collective

     55:   Input Parameter:
     56: . comm - The communicator for the `DMMOAB` object

     58:   Output Parameter:
     59: . dmb - The `DMMOAB` object

     61:   Level: beginner

     63: .seealso: `DMMOAB`, `DMMoabCreateMoab()`
     64: @*/
     65: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
     66: {
     67:   PetscFunctionBegin;
     68:   PetscAssertPointer(dmb, 2);
     69:   PetscCall(DMCreate(comm, dmb));
     70:   PetscCall(DMSetType(*dmb, DMMOAB));
     71:   PetscFunctionReturn(PETSC_SUCCESS);
     72: }

     74: /*@C
     75:   DMMoabCreateMoab - Creates a `DMMOAB` object, optionally from an instance and other data

     77:   Collective

     79:   Input Parameters:
     80: + comm     - The communicator for the `DMMOAB` object
     81: . mbiface  - (ptr to) the `DMMOAB` Instance; if passed in `NULL`, MOAB instance is created inside PETSc, and destroyed
     82:          along with the `DMMOAB`
     83: . ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
     84: - range    - If non-`NULL`, contains range of entities to which DOFs will be assigned

     86:   Output Parameter:
     87: . dmb - The `DMMOAB` object

     89:   Level: intermediate

     91: .seealso: `DMMOAB`, `DMMoabCreate()`
     92: @*/
     93: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
     94: {
     95:   moab::ErrorCode merr;
     96:   DM              dmmb;
     97:   DM_Moab        *dmmoab;

     99:   PetscFunctionBegin;
    100:   PetscAssertPointer(dmb, 6);

    102:   PetscCall(DMMoabCreate(comm, &dmmb));
    103:   dmmoab = (DM_Moab *)dmmb->data;

    105:   if (!mbiface) {
    106:     dmmoab->mbiface          = new moab::Core();
    107:     dmmoab->icreatedinstance = PETSC_TRUE;
    108:   } else {
    109:     dmmoab->mbiface          = mbiface;
    110:     dmmoab->icreatedinstance = PETSC_FALSE;
    111:   }

    113:   /* by default the fileset = root set. This set stores the hierarchy of entities belonging to current DM */
    114:   dmmoab->fileset     = 0;
    115:   dmmoab->hlevel      = 0;
    116:   dmmoab->nghostrings = 0;

    118: #ifdef MOAB_HAVE_MPI
    119:   moab::EntityHandle partnset;

    121:   /* Create root sets for each mesh.  Then pass these
    122:       to the load_file functions to be populated. */
    123:   merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, partnset);
    124:   MBERR("Creating partition set failed", merr);

    126:   /* Create the parallel communicator object with the partition handle associated with MOAB */
    127:   dmmoab->pcomm = moab::ParallelComm::get_pcomm(dmmoab->mbiface, partnset, &comm);
    128: #endif

    130:   /* do the remaining initializations for DMMoab */
    131:   dmmoab->bs        = 1;
    132:   dmmoab->numFields = 1;
    133:   PetscCall(PetscMalloc(dmmoab->numFields * sizeof(char *), &dmmoab->fieldNames));
    134:   PetscCall(PetscStrallocpy("DEFAULT", (char **)&dmmoab->fieldNames[0]));
    135:   dmmoab->rw_dbglevel            = 0;
    136:   dmmoab->partition_by_rank      = PETSC_FALSE;
    137:   dmmoab->extra_read_options[0]  = '\0';
    138:   dmmoab->extra_write_options[0] = '\0';
    139:   dmmoab->read_mode              = READ_PART;
    140:   dmmoab->write_mode             = WRITE_PART;

    142:   /* set global ID tag handle */
    143:   if (ltog_tag && *ltog_tag) {
    144:     PetscCall(DMMoabSetLocalToGlobalTag(dmmb, *ltog_tag));
    145:   } else {
    146:     merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
    147:     MBERRNM(merr);
    148:     if (ltog_tag) *ltog_tag = dmmoab->ltog_tag;
    149:   }

    151:   merr = dmmoab->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dmmoab->material_tag);
    152:   MBERRNM(merr);

    154:   /* set the local range of entities (vertices) of interest */
    155:   if (range) PetscCall(DMMoabSetLocalVertices(dmmb, range));
    156:   *dmb = dmmb;
    157:   PetscFunctionReturn(PETSC_SUCCESS);
    158: }

    160: #ifdef MOAB_HAVE_MPI

    162: /*@C
    163:   DMMoabGetParallelComm - Get the ParallelComm used with this `DMMOAB`

    165:   Collective

    167:   Input Parameter:
    168: . dm    - The `DMMOAB` object being set

    170:   Output Parameter:
    171: . pcomm - The ParallelComm for the `DMMOAB`

    173:   Level: beginner

    175: .seealso: `DMMOAB`, `DMMoabSetInterface()`
    176: @*/
    177: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
    178: {
    179:   PetscFunctionBegin;
    181:   *pcomm = ((DM_Moab *)dm->data)->pcomm;
    182:   PetscFunctionReturn(PETSC_SUCCESS);
    183: }

    185: #endif /* MOAB_HAVE_MPI */

    187: /*@C
    188:   DMMoabSetInterface - Set the MOAB instance used with this `DMMOAB`

    190:   Collective

    192:   Input Parameters:
    193: + dm      - The `DMMOAB` object being set
    194: - mbiface - The MOAB instance being set on this `DMMOAB`

    196:   Level: beginner

    198: .seealso: `DMMOAB`, `DMMoabGetInterface()`
    199: @*/
    200: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
    201: {
    202:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

    204:   PetscFunctionBegin;
    206:   PetscAssertPointer(mbiface, 2);
    207: #ifdef MOAB_HAVE_MPI
    208:   dmmoab->pcomm = NULL;
    209: #endif
    210:   dmmoab->mbiface          = mbiface;
    211:   dmmoab->icreatedinstance = PETSC_FALSE;
    212:   PetscFunctionReturn(PETSC_SUCCESS);
    213: }

    215: /*@C
    216:   DMMoabGetInterface - Get the MOAB instance used with this `DMMOAB`

    218:   Collective

    220:   Input Parameter:
    221: . dm - The `DMMOAB` object being set

    223:   Output Parameter:
    224: . mbiface - The MOAB instance set on this `DMMOAB`

    226:   Level: beginner

    228: .seealso: `DMMOAB`, `DMMoabSetInterface()`
    229: @*/
    230: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
    231: {
    232:   static PetscBool cite = PETSC_FALSE;

    234:   PetscFunctionBegin;
    236:   PetscCall(PetscCitationsRegister("@techreport{tautges_moab:_2004,\n  type = {{SAND2004-1592}},\n  title = {{MOAB:} A Mesh-Oriented Database},  institution = {Sandia National Laboratories},\n  author = {Tautges, T. J. and Meyers, R. and Merkley, "
    237:                                    "K. and Stimpson, C. and Ernst, C.},\n  year = {2004},  note = {Report}\n}\n",
    238:                                    &cite));
    239:   *mbiface = ((DM_Moab *)dm->data)->mbiface;
    240:   PetscFunctionReturn(PETSC_SUCCESS);
    241: }

    243: /*@C
    244:   DMMoabSetLocalVertices - Set the entities having DOFs on this `DMMOAB`

    246:   Collective

    248:   Input Parameters:
    249: + dm    - The `DMMOAB` object being set
    250: - range - The entities treated by this `DMMOAB`

    252:   Level: beginner

    254: .seealso: `DMMOAB`, `DMMoabGetAllVertices()`
    255: @*/
    256: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
    257: {
    258:   moab::Range tmpvtxs;
    259:   DM_Moab    *dmmoab = (DM_Moab *)dm->data;

    261:   PetscFunctionBegin;
    263:   dmmoab->vlocal->clear();
    264:   dmmoab->vowned->clear();

    266:   dmmoab->vlocal->insert(range->begin(), range->end());

    268: #ifdef MOAB_HAVE_MPI
    269:   moab::ErrorCode merr;
    270:   /* filter based on parallel status */
    271:   merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
    272:   MBERRNM(merr);

    274:   /* filter all the non-owned and shared entities out of the list */
    275:   tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
    276:   merr    = dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
    277:   MBERRNM(merr);
    278:   tmpvtxs         = moab::subtract(tmpvtxs, *dmmoab->vghost);
    279:   *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
    280: #else
    281:   *dmmoab->vowned = *dmmoab->vlocal;
    282: #endif

    284:   /* compute and cache the sizes of local and ghosted entities */
    285:   dmmoab->nloc   = dmmoab->vowned->size();
    286:   dmmoab->nghost = dmmoab->vghost->size();
    287: #ifdef MOAB_HAVE_MPI
    288:   PetscCallMPI(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
    289: #else
    290:   dmmoab->n = dmmoab->nloc;
    291: #endif
    292:   PetscFunctionReturn(PETSC_SUCCESS);
    293: }

    295: /*@C
    296:   DMMoabGetAllVertices - Get the entities having DOFs on this `DMMOAB`

    298:   Collective

    300:   Input Parameter:
    301: . dm - The `DMMOAB` object being set

    303:   Output Parameter:
    304: . local - The local vertex entities in this `DMMOAB` = (owned+ghosted)

    306:   Level: beginner

    308: .seealso: `DMMOAB`, `DMMoabGetLocalVertices()`
    309: @*/
    310: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
    311: {
    312:   PetscFunctionBegin;
    314:   if (local) *local = *((DM_Moab *)dm->data)->vlocal;
    315:   PetscFunctionReturn(PETSC_SUCCESS);
    316: }

    318: /*@C
    319:   DMMoabGetLocalVertices - Get the entities having DOFs on this `DMMOAB`

    321:   Collective

    323:   Input Parameter:
    324: . dm - The `DMMOAB` object being set

    326:   Output Parameters:
    327: + owned - The owned vertex entities in this `DMMOAB`
    328: - ghost - The ghosted entities (non-owned) stored locally in this partition

    330:   Level: beginner

    332: .seealso: `DMMOAB`, `DMMoabGetAllVertices()`
    333: @*/
    334: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
    335: {
    336:   PetscFunctionBegin;
    338:   if (owned) *owned = ((DM_Moab *)dm->data)->vowned;
    339:   if (ghost) *ghost = ((DM_Moab *)dm->data)->vghost;
    340:   PetscFunctionReturn(PETSC_SUCCESS);
    341: }

    343: /*@C
    344:   DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned

    346:   Collective

    348:   Input Parameter:
    349: . dm - The `DMMOAB` object being set

    351:   Output Parameter:
    352: . range - The entities owned locally

    354:   Level: beginner

    356: .seealso: `DMMOAB`, `DMMoabSetLocalElements()`
    357: @*/
    358: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
    359: {
    360:   PetscFunctionBegin;
    362:   if (range) *range = ((DM_Moab *)dm->data)->elocal;
    363:   PetscFunctionReturn(PETSC_SUCCESS);
    364: }

    366: /*@C
    367:   DMMoabSetLocalElements - Set the entities having DOFs on this `DMMOAB`

    369:   Collective

    371:   Input Parameters:
    372: + dm    - The `DMMOAB` object being set
    373: - range - The entities treated by this `DMMOAB`

    375:   Level: beginner

    377: .seealso: `DMMOAB`, `DMMoabGetLocalElements()`
    378: @*/
    379: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
    380: {
    381:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

    383:   PetscFunctionBegin;
    385:   dmmoab->elocal->clear();
    386:   dmmoab->eghost->clear();
    387:   dmmoab->elocal->insert(range->begin(), range->end());
    388: #ifdef MOAB_HAVE_MPI
    389:   moab::ErrorCode merr;
    390:   merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
    391:   MBERRNM(merr);
    392:   *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
    393: #endif
    394:   dmmoab->neleloc   = dmmoab->elocal->size();
    395:   dmmoab->neleghost = dmmoab->eghost->size();
    396: #ifdef MOAB_HAVE_MPI
    397:   PetscCallMPI(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
    398:   PetscCall(PetscInfo(dm, "Created %" PetscInt_FMT " local and %" PetscInt_FMT " global elements.\n", dmmoab->neleloc, dmmoab->nele));
    399: #else
    400:   dmmoab->nele = dmmoab->neleloc;
    401: #endif
    402:   PetscFunctionReturn(PETSC_SUCCESS);
    403: }

    405: /*@C
    406:   DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering

    408:   Collective

    410:   Input Parameters:
    411: + dm      - The `DMMOAB` object being set
    412: - ltogtag - The `DMMOAB` tag used for local to global ids

    414:   Level: beginner

    416: .seealso: `DMMOAB`, `DMMoabGetLocalToGlobalTag()`
    417: @*/
    418: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
    419: {
    420:   PetscFunctionBegin;
    422:   ((DM_Moab *)dm->data)->ltog_tag = ltogtag;
    423:   PetscFunctionReturn(PETSC_SUCCESS);
    424: }

    426: /*@C
    427:   DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering

    429:   Collective

    431:   Input Parameter:
    432: . dm - The `DMMOAB` object being set

    434:   Output Parameter:
    435: . ltog_tag - The MOAB tag used for local to global ids

    437:   Level: beginner

    439: .seealso: `DMMOAB`, `DMMoabSetLocalToGlobalTag()`
    440: @*/
    441: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
    442: {
    443:   PetscFunctionBegin;
    445:   *ltog_tag = ((DM_Moab *)dm->data)->ltog_tag;
    446:   PetscFunctionReturn(PETSC_SUCCESS);
    447: }

    449: /*@C
    450:   DMMoabSetBlockSize - Set the block size used with this `DMMOAB`

    452:   Collective

    454:   Input Parameters:
    455: + dm - The `DMMOAB` object being set
    456: - bs - The block size used with this `DMMOAB`

    458:   Level: beginner

    460: .seealso: `DMMOAB`, `DMMoabGetBlockSize()`
    461: @*/
    462: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
    463: {
    464:   PetscFunctionBegin;
    466:   ((DM_Moab *)dm->data)->bs = bs;
    467:   PetscFunctionReturn(PETSC_SUCCESS);
    468: }

    470: /*@C
    471:   DMMoabGetBlockSize - Get the block size used with this `DMMOAB`

    473:   Collective

    475:   Input Parameter:
    476: . dm - The `DMMOAB` object being set

    478:   Output Parameter:
    479: . bs - The block size used with this `DMMOAB`

    481:   Level: beginner

    483: .seealso: `DMMOAB`, `DMMoabSetBlockSize()`
    484: @*/
    485: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
    486: {
    487:   PetscFunctionBegin;
    489:   *bs = ((DM_Moab *)dm->data)->bs;
    490:   PetscFunctionReturn(PETSC_SUCCESS);
    491: }

    493: /*@C
    494:   DMMoabGetSize - Get the global vertex size used with this `DMMOAB`

    496:   Collective

    498:   Input Parameter:
    499: . dm - The `DMMOAB` object being set

    501:   Output Parameters:
    502: + neg - The number of global elements in the `DMMOAB` instance
    503: - nvg - The number of global vertices in the `DMMOAB` instance

    505:   Level: beginner

    507: .seealso: `DMMOAB`, `DMMoabGetLocalSize()`
    508: @*/
    509: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
    510: {
    511:   PetscFunctionBegin;
    513:   if (neg) *neg = ((DM_Moab *)dm->data)->nele;
    514:   if (nvg) *nvg = ((DM_Moab *)dm->data)->n;
    515:   PetscFunctionReturn(PETSC_SUCCESS);
    516: }

    518: /*@C
    519:   DMMoabGetLocalSize - Get the local and ghosted vertex size used with this `DMMOAB`

    521:   Collective

    523:   Input Parameter:
    524: . dm - The `DMMOAB` object being set

    526:   Output Parameters:
    527: + nel - The number of owned elements in this processor
    528: . neg - The number of ghosted elements in this processor
    529: . nvl - The number of owned vertices in this processor
    530: - nvg - The number of ghosted vertices in this processor

    532:   Level: beginner

    534: .seealso: `DMMOAB`, `DMMoabGetSize()`
    535: @*/
    536: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
    537: {
    538:   PetscFunctionBegin;
    540:   if (nel) *nel = ((DM_Moab *)dm->data)->neleloc;
    541:   if (neg) *neg = ((DM_Moab *)dm->data)->neleghost;
    542:   if (nvl) *nvl = ((DM_Moab *)dm->data)->nloc;
    543:   if (nvg) *nvg = ((DM_Moab *)dm->data)->nghost;
    544:   PetscFunctionReturn(PETSC_SUCCESS);
    545: }

    547: /*@C
    548:   DMMoabGetOffset - Get the local offset for the global vector

    550:   Collective

    552:   Input Parameter:
    553: . dm - The `DMMOAB` object being set

    555:   Output Parameter:
    556: . offset - The local offset for the global vector

    558:   Level: beginner

    560: .seealso: `DMMOAB`, `DMMoabGetDimension()`
    561: @*/
    562: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
    563: {
    564:   PetscFunctionBegin;
    566:   *offset = ((DM_Moab *)dm->data)->vstart;
    567:   PetscFunctionReturn(PETSC_SUCCESS);
    568: }

    570: /*@C
    571:   DMMoabGetDimension - Get the dimension of the `DM` Mesh

    573:   Collective

    575:   Input Parameter:
    576: . dm - The `DMMOAB` object

    578:   Output Parameter:
    579: . dim - The dimension of `DM`

    581:   Level: beginner

    583: .seealso: `DMMOAB`, `DMMoabGetOffset()`
    584: @*/
    585: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
    586: {
    587:   PetscFunctionBegin;
    589:   *dim = ((DM_Moab *)dm->data)->dim;
    590:   PetscFunctionReturn(PETSC_SUCCESS);
    591: }

    593: /*@C
    594:   DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
    595:   generated through uniform refinement.

    597:   Collective

    599:   Input Parameter:
    600: . dm - The `DMMOAB` object being set

    602:   Output Parameter:
    603: . nlevel - The current mesh hierarchy level

    605:   Level: beginner

    607: .seealso: `DMMOAB`, `DMMoabGetMaterialBlock()`
    608: @*/
    609: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
    610: {
    611:   PetscFunctionBegin;
    613:   if (nlevel) *nlevel = ((DM_Moab *)dm->data)->hlevel;
    614:   PetscFunctionReturn(PETSC_SUCCESS);
    615: }

    617: /*@C
    618:   DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh

    620:   Collective

    622:   Input Parameters:
    623: + dm      - The `DMMOAB` object
    624: - ehandle - The element entity handle

    626:   Output Parameter:
    627: . mat - The material ID for the current entity

    629:   Level: beginner

    631: .seealso: `DMMOAB`, `DMMoabGetHierarchyLevel()`
    632: @*/
    633: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
    634: {
    635:   DM_Moab *dmmoab;

    637:   PetscFunctionBegin;
    639:   if (*mat) {
    640:     dmmoab = (DM_Moab *)dm->data;
    641:     *mat   = dmmoab->materials[dmmoab->elocal->index(ehandle)];
    642:   }
    643:   PetscFunctionReturn(PETSC_SUCCESS);
    644: }

    646: /*@C
    647:   DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities

    649:   Collective

    651:   Input Parameters:
    652: + dm    - The `DMMOAB` object
    653: . nconn - Number of entities whose coordinates are needed
    654: - conn  - The vertex entity handles

    656:   Output Parameter:
    657: . vpos - The coordinates of the requested vertex entities

    659:   Level: beginner

    661: .seealso: `DMMOAB`, `DMMoabGetVertexConnectivity()`
    662: @*/
    663: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
    664: {
    665:   DM_Moab        *dmmoab;
    666:   moab::ErrorCode merr;

    668:   PetscFunctionBegin;
    670:   PetscAssertPointer(conn, 3);
    671:   PetscAssertPointer(vpos, 4);
    672:   dmmoab = (DM_Moab *)dm->data;

    674:   /* Get connectivity information in MOAB canonical ordering */
    675:   if (dmmoab->hlevel) {
    676:     merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle *>(conn), nconn, dmmoab->hlevel, vpos);
    677:     MBERRNM(merr);
    678:   } else {
    679:     merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);
    680:     MBERRNM(merr);
    681:   }
    682:   PetscFunctionReturn(PETSC_SUCCESS);
    683: }

    685: /*@C
    686:   DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity

    688:   Collective

    690:   Input Parameters:
    691: + dm      - The `DMMOAB` object
    692: - vhandle - Vertex entity handle

    694:   Output Parameters:
    695: + nconn - Number of entities whose coordinates are needed
    696: - conn  - The vertex entity handles

    698:   Level: beginner

    700: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabRestoreVertexConnectivity()`
    701: @*/
    702: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt *nconn, moab::EntityHandle **conn)
    703: {
    704:   DM_Moab                        *dmmoab;
    705:   std::vector<moab::EntityHandle> adj_entities, connect;
    706:   moab::ErrorCode                 merr;

    708:   PetscFunctionBegin;
    710:   PetscAssertPointer(conn, 4);
    711:   dmmoab = (DM_Moab *)dm->data;

    713:   /* Get connectivity information in MOAB canonical ordering */
    714:   merr = dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION);
    715:   MBERRNM(merr);
    716:   merr = dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect);
    717:   MBERRNM(merr);

    719:   if (conn) {
    720:     PetscCall(PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn));
    721:     PetscCall(PetscArraycpy(*conn, &connect[0], connect.size()));
    722:   }
    723:   if (nconn) *nconn = connect.size();
    724:   PetscFunctionReturn(PETSC_SUCCESS);
    725: }

    727: /*@C
    728:   DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity

    730:   Collective

    732:   Input Parameters:
    733: + dm      - The `DMMOAB` object
    734: . ehandle - Vertex entity handle
    735: . nconn   - Number of entities whose coordinates are needed
    736: - conn    - The vertex entity handles

    738:   Level: beginner

    740: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`
    741: @*/
    742: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, moab::EntityHandle **conn)
    743: {
    744:   PetscFunctionBegin;
    746:   PetscAssertPointer(conn, 4);

    748:   if (conn) PetscCall(PetscFree(*conn));
    749:   if (nconn) *nconn = 0;
    750:   PetscFunctionReturn(PETSC_SUCCESS);
    751: }

    753: /*@C
    754:   DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity

    756:   Collective

    758:   Input Parameters:
    759: + dm      - The `DMMOAB` object
    760: - ehandle - Vertex entity handle

    762:   Output Parameters:
    763: + nconn - Number of entities whose coordinates are needed
    764: - conn  - The vertex entity handles

    766:   Level: beginner

    768: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`, `DMMoabRestoreVertexConnectivity()`
    769: @*/
    770: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, const moab::EntityHandle **conn)
    771: {
    772:   DM_Moab                        *dmmoab;
    773:   const moab::EntityHandle       *connect;
    774:   std::vector<moab::EntityHandle> vconn;
    775:   moab::ErrorCode                 merr;
    776:   PetscInt                        nnodes;

    778:   PetscFunctionBegin;
    780:   PetscAssertPointer(conn, 4);
    781:   dmmoab = (DM_Moab *)dm->data;

    783:   /* Get connectivity information in MOAB canonical ordering */
    784:   merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes);
    785:   MBERRNM(merr);
    786:   if (conn) *conn = connect;
    787:   if (nconn) *nconn = nnodes;
    788:   PetscFunctionReturn(PETSC_SUCCESS);
    789: }

    791: /*@C
    792:   DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)

    794:   Collective

    796:   Input Parameters:
    797: + dm  - The `DMMOAB` object
    798: - ent - Entity handle

    800:   Output Parameter:
    801: . ent_on_boundary - `PETSC_TRUE` if entity on boundary; `PETSC_FALSE` otherwise

    803:   Level: beginner

    805: .seealso: `DMMOAB`, `DMMoabCheckBoundaryVertices()`
    806: @*/
    807: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool *ent_on_boundary)
    808: {
    809:   moab::EntityType etype;
    810:   DM_Moab         *dmmoab;
    811:   PetscInt         edim;

    813:   PetscFunctionBegin;
    815:   PetscAssertPointer(ent_on_boundary, 3);
    816:   dmmoab = (DM_Moab *)dm->data;

    818:   /* get the entity type and handle accordingly */
    819:   etype = dmmoab->mbiface->type_from_handle(ent);
    820:   PetscCheck(etype < moab::MBPOLYHEDRON, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Entity type on the boundary skin is invalid. EntityType = %" PetscInt_FMT, etype);

    822:   /* get the entity dimension */
    823:   edim = dmmoab->mbiface->dimension_from_handle(ent);

    825:   *ent_on_boundary = PETSC_FALSE;
    826:   if (etype == moab::MBVERTEX && edim == 0) {
    827:     *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
    828:   } else {
    829:     if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
    830:       if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
    831:     } else { /* next check the lower-dimensional faces */
    832:       if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
    833:     }
    834:   }
    835:   PetscFunctionReturn(PETSC_SUCCESS);
    836: }

    838: /*@C
    839:   DMMoabCheckBoundaryVertices - Check whether a given entity is on the boundary (vertex, edge, face, element)

    841:   Input Parameters:
    842: + dm    - The `DMMOAB` object
    843: . nconn - Number of handles
    844: - cnt   - Array of entity handles

    846:   Output Parameter:
    847: . isbdvtx - Array of boundary markers - `PETSC_TRUE` if entity on boundary; `PETSC_FALSE` otherwise

    849:   Level: beginner

    851: .seealso: `DMMOAB`, `DMMoabIsEntityOnBoundary()`
    852: @*/
    853: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool *isbdvtx)
    854: {
    855:   DM_Moab *dmmoab;
    856:   PetscInt i;

    858:   PetscFunctionBegin;
    860:   PetscAssertPointer(cnt, 3);
    861:   PetscAssertPointer(isbdvtx, 4);
    862:   dmmoab = (DM_Moab *)dm->data;

    864:   for (i = 0; i < nconn; ++i) isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
    865:   PetscFunctionReturn(PETSC_SUCCESS);
    866: }

    868: /*@C
    869:   DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary

    871:   Input Parameter:
    872: . dm - The `DMMOAB` object

    874:   Output Parameters:
    875: + bdvtx   - Boundary vertices
    876: . bdelems - Boundary elements
    877: - bdfaces - Boundary faces

    879:   Level: beginner

    881: .seealso: `DMMOAB`, `DMMoabCheckBoundaryVertices()`, `DMMoabIsEntityOnBoundary()`
    882: @*/
    883: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range **bdelems, const moab::Range **bdfaces)
    884: {
    885:   DM_Moab *dmmoab;

    887:   PetscFunctionBegin;
    889:   dmmoab = (DM_Moab *)dm->data;

    891:   if (bdvtx) *bdvtx = dmmoab->bndyvtx;
    892:   if (bdfaces) *bdfaces = dmmoab->bndyfaces;
    893:   if (bdelems) *bdfaces = dmmoab->bndyelems;
    894:   PetscFunctionReturn(PETSC_SUCCESS);
    895: }

    897: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
    898: {
    899:   PetscInt        i;
    900:   moab::ErrorCode merr;
    901:   DM_Moab        *dmmoab = (DM_Moab *)dm->data;

    903:   PetscFunctionBegin;

    906:   dmmoab->refct--;
    907:   if (!dmmoab->refct) {
    908:     delete dmmoab->vlocal;
    909:     delete dmmoab->vowned;
    910:     delete dmmoab->vghost;
    911:     delete dmmoab->elocal;
    912:     delete dmmoab->eghost;
    913:     delete dmmoab->bndyvtx;
    914:     delete dmmoab->bndyfaces;
    915:     delete dmmoab->bndyelems;

    917:     PetscCall(PetscFree(dmmoab->gsindices));
    918:     PetscCall(PetscFree2(dmmoab->gidmap, dmmoab->lidmap));
    919:     PetscCall(PetscFree(dmmoab->dfill));
    920:     PetscCall(PetscFree(dmmoab->ofill));
    921:     PetscCall(PetscFree(dmmoab->materials));
    922:     if (dmmoab->fieldNames) {
    923:       for (i = 0; i < dmmoab->numFields; i++) PetscCall(PetscFree(dmmoab->fieldNames[i]));
    924:       PetscCall(PetscFree(dmmoab->fieldNames));
    925:     }

    927:     if (dmmoab->nhlevels) {
    928:       PetscCall(PetscFree(dmmoab->hsets));
    929:       dmmoab->nhlevels = 0;
    930:       if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
    931:       dmmoab->hierarchy = NULL;
    932:     }

    934:     if (dmmoab->icreatedinstance) {
    935:       delete dmmoab->pcomm;
    936:       merr = dmmoab->mbiface->delete_mesh();
    937:       MBERRNM(merr);
    938:       delete dmmoab->mbiface;
    939:     }
    940:     dmmoab->mbiface = NULL;
    941: #ifdef MOAB_HAVE_MPI
    942:     dmmoab->pcomm = NULL;
    943: #endif
    944:     PetscCall(VecScatterDestroy(&dmmoab->ltog_sendrecv));
    945:     PetscCall(ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map));
    946:     PetscCall(PetscFree(dm->data));
    947:   }
    948:   PetscFunctionReturn(PETSC_SUCCESS);
    949: }

    951: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(DM dm, PetscOptionItems *PetscOptionsObject)
    952: {
    953:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

    955:   PetscFunctionBegin;
    956:   PetscOptionsHeadBegin(PetscOptionsObject, "DMMoab Options");
    957:   PetscCall(PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL, 0));
    958:   PetscCall(PetscOptionsBool("-dm_moab_partiton_by_rank", "Use partition by rank when reading MOAB meshes from file", "DMView", dmmoab->partition_by_rank, &dmmoab->partition_by_rank, NULL));
    959:   /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
    960:   PetscCall(PetscOptionsString("-dm_moab_read_opts", "Extra options to enable MOAB reader to load DM from file", "DMView", dmmoab->extra_read_options, dmmoab->extra_read_options, sizeof(dmmoab->extra_read_options), NULL));
    961:   PetscCall(PetscOptionsString("-dm_moab_write_opts", "Extra options to enable MOAB writer to serialize DM to file", "DMView", dmmoab->extra_write_options, dmmoab->extra_write_options, sizeof(dmmoab->extra_write_options), NULL));
    962:   PetscCall(PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum *)&dmmoab->read_mode, NULL));
    963:   PetscCall(PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum *)&dmmoab->write_mode, NULL));
    964:   PetscOptionsHeadEnd();
    965:   PetscFunctionReturn(PETSC_SUCCESS);
    966: }

    968: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
    969: {
    970:   moab::ErrorCode       merr;
    971:   Vec                   local, global;
    972:   IS                    from, to;
    973:   moab::Range::iterator iter;
    974:   PetscInt              i, j, f, bs, vent, totsize, *lgmap;
    975:   DM_Moab              *dmmoab = (DM_Moab *)dm->data;
    976:   moab::Range           adjs;

    978:   PetscFunctionBegin;
    980:   /* Get the local and shared vertices and cache it */
    981:   PetscCheck(dmmoab->mbiface != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB Interface before calling SetUp.");
    982: #ifdef MOAB_HAVE_MPI
    983:   PetscCheck(dmmoab->pcomm != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB ParallelComm object before calling SetUp.");
    984: #endif

    986:   /* Get the entities recursively in the current part of the mesh, if user did not set the local vertices explicitly */
    987:   if (dmmoab->vlocal->empty()) {
    988:     //merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,*dmmoab->vlocal,true);MBERRNM(merr);
    989:     merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, 0, *dmmoab->vlocal, false);
    990:     MBERRNM(merr);

    992: #ifdef MOAB_HAVE_MPI
    993:     /* filter based on parallel status */
    994:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
    995:     MBERRNM(merr);

    997:     /* filter all the non-owned and shared entities out of the list */
    998:     // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
    999:     adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
    1000:     merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
    1001:     MBERRNM(merr);
    1002:     adjs            = moab::subtract(adjs, *dmmoab->vghost);
    1003:     *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
    1004: #else
    1005:     *dmmoab->vowned = *dmmoab->vlocal;
    1006: #endif

    1008:     /* compute and cache the sizes of local and ghosted entities */
    1009:     dmmoab->nloc   = dmmoab->vowned->size();
    1010:     dmmoab->nghost = dmmoab->vghost->size();

    1012: #ifdef MOAB_HAVE_MPI
    1013:     PetscCallMPI(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
    1014:     PetscCall(PetscInfo(NULL, "Filset ID: %lu, Vertices: local - %zu, owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->fileset, dmmoab->vlocal->size(), dmmoab->nloc, dmmoab->nghost));
    1015: #else
    1016:     dmmoab->n = dmmoab->nloc;
    1017: #endif
    1018:   }

    1020:   {
    1021:     /* get the information about the local elements in the mesh */
    1022:     dmmoab->eghost->clear();

    1024:     /* first decipher the leading dimension */
    1025:     for (i = 3; i > 0; i--) {
    1026:       dmmoab->elocal->clear();
    1027:       merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false);
    1028:       MBERRNM(merr);

    1030:       /* store the current mesh dimension */
    1031:       if (dmmoab->elocal->size()) {
    1032:         dmmoab->dim = i;
    1033:         break;
    1034:       }
    1035:     }

    1037:     PetscCall(DMSetDimension(dm, dmmoab->dim));

    1039: #ifdef MOAB_HAVE_MPI
    1040:     /* filter the ghosted and owned element list */
    1041:     *dmmoab->eghost = *dmmoab->elocal;
    1042:     merr            = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
    1043:     MBERRNM(merr);
    1044:     *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
    1045: #endif

    1047:     dmmoab->neleloc   = dmmoab->elocal->size();
    1048:     dmmoab->neleghost = dmmoab->eghost->size();

    1050: #ifdef MOAB_HAVE_MPI
    1051:     PetscCallMPI(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
    1052:     PetscCall(PetscInfo(NULL, "%d-dim elements: owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost));
    1053: #else
    1054:     dmmoab->nele = dmmoab->neleloc;
    1055: #endif
    1056:   }

    1058:   bs = dmmoab->bs;
    1059:   if (!dmmoab->ltog_tag) {
    1060:     /* Get the global ID tag. The global ID tag is applied to each
    1061:        vertex. It acts as an global identifier which MOAB uses to
    1062:        assemble the individual pieces of the mesh */
    1063:     merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
    1064:     MBERRNM(merr);
    1065:   }

    1067:   totsize = dmmoab->vlocal->size();
    1068:   PetscCheck(totsize == dmmoab->nloc + dmmoab->nghost, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Mismatch between local and owned+ghost vertices. %" PetscInt_FMT " != %" PetscInt_FMT ".", totsize, dmmoab->nloc + dmmoab->nghost);
    1069:   PetscCall(PetscCalloc1(totsize, &dmmoab->gsindices));
    1070:   {
    1071:     /* first get the local indices */
    1072:     merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]);
    1073:     MBERRNM(merr);
    1074:     if (dmmoab->nghost) { /* next get the ghosted indices */
    1075:       merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]);
    1076:       MBERRNM(merr);
    1077:     }

    1079:     /* find out the local and global minima of GLOBAL_ID */
    1080:     dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
    1081:     for (i = 0; i < totsize; ++i) {
    1082:       if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
    1083:       if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
    1084:     }

    1086:     PetscCallMPI(MPIU_Allreduce(&dmmoab->lminmax[0], &dmmoab->gminmax[0], 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm));
    1087:     PetscCallMPI(MPIU_Allreduce(&dmmoab->lminmax[1], &dmmoab->gminmax[1], 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm));

    1089:     /* set the GID map */
    1090:     for (i = 0; i < totsize; ++i) { dmmoab->gsindices[i] -= dmmoab->gminmax[0]; /* zero based index needed for IS */ }
    1091:     dmmoab->lminmax[0] -= dmmoab->gminmax[0];
    1092:     dmmoab->lminmax[1] -= dmmoab->gminmax[0];

    1094:     PetscCall(PetscInfo(NULL, "GLOBAL_ID: Local [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "], Global [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "]\n", dmmoab->lminmax[0], dmmoab->lminmax[1], dmmoab->gminmax[0], dmmoab->gminmax[1]));
    1095:   }
    1096:   PetscCheck(dmmoab->bs == dmmoab->numFields || dmmoab->bs == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mismatch between block size and number of component fields. %" PetscInt_FMT " != 1 OR %" PetscInt_FMT " != %" PetscInt_FMT ".", dmmoab->bs, dmmoab->bs,
    1097:              dmmoab->numFields);

    1099:   {
    1100:     dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
    1101:     dmmoab->seqend   = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
    1102:     PetscCall(PetscInfo(NULL, "SEQUENCE: Local [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "]\n", dmmoab->seqstart, dmmoab->seqend));

    1104:     PetscCall(PetscMalloc2(dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->gidmap, dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->lidmap));
    1105:     PetscCall(PetscMalloc1(totsize * dmmoab->numFields, &lgmap));

    1107:     i = j = 0;
    1108:     /* set the owned vertex data first */
    1109:     for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
    1110:       vent                 = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
    1111:       dmmoab->gidmap[vent] = dmmoab->gsindices[i];
    1112:       dmmoab->lidmap[vent] = i;
    1113:       for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
    1114:     }
    1115:     /* next arrange all the ghosted data information */
    1116:     for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
    1117:       vent                 = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
    1118:       dmmoab->gidmap[vent] = dmmoab->gsindices[i];
    1119:       dmmoab->lidmap[vent] = i;
    1120:       for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
    1121:     }

    1123:     /* We need to create the Global to Local Vector Scatter Contexts
    1124:        1) First create a local and global vector
    1125:        2) Create a local and global IS
    1126:        3) Create VecScatter and LtoGMapping objects
    1127:        4) Cleanup the IS and Vec objects
    1128:     */
    1129:     PetscCall(DMCreateGlobalVector(dm, &global));
    1130:     PetscCall(DMCreateLocalVector(dm, &local));

    1132:     PetscCall(VecGetOwnershipRange(global, &dmmoab->vstart, &dmmoab->vend));

    1134:     /* global to local must retrieve ghost points */
    1135:     PetscCall(ISCreateStride(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, dmmoab->vstart, 1, &from));
    1136:     PetscCall(ISSetBlockSize(from, bs));

    1138:     PetscCall(ISCreateGeneral(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, &lgmap[0], PETSC_COPY_VALUES, &to));
    1139:     PetscCall(ISSetBlockSize(to, bs));

    1141:     if (!dmmoab->ltog_map) {
    1142:       /* create to the local to global mapping for vectors in order to use VecSetValuesLocal */
    1143:       PetscCall(ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm, dmmoab->bs, totsize * dmmoab->numFields, lgmap, PETSC_COPY_VALUES, &dmmoab->ltog_map));
    1144:     }

    1146:     /* now create the scatter object from local to global vector */
    1147:     PetscCall(VecScatterCreate(local, from, global, to, &dmmoab->ltog_sendrecv));

    1149:     /* clean up IS, Vec */
    1150:     PetscCall(PetscFree(lgmap));
    1151:     PetscCall(ISDestroy(&from));
    1152:     PetscCall(ISDestroy(&to));
    1153:     PetscCall(VecDestroy(&local));
    1154:     PetscCall(VecDestroy(&global));
    1155:   }

    1157:   dmmoab->bndyvtx   = new moab::Range();
    1158:   dmmoab->bndyfaces = new moab::Range();
    1159:   dmmoab->bndyelems = new moab::Range();
    1160:   /* skin the boundary and store nodes */
    1161:   if (!dmmoab->hlevel) {
    1162:     /* get the skin vertices of boundary faces for the current partition and then filter
    1163:        the local, boundary faces, vertices and elements alone via PSTATUS flags;
    1164:        this should not give us any ghosted boundary, but if user needs such a functionality
    1165:        it would be easy to add it based on the find_skin query below */
    1166:     moab::Skinner skinner(dmmoab->mbiface);

    1168:     /* get the entities on the skin - only the faces */
    1169:     merr = skinner.find_skin(dmmoab->fileset, *dmmoab->elocal, false, *dmmoab->bndyfaces, NULL, true, true, false);
    1170:     MBERRNM(merr); // 'false' param indicates we want faces back, not vertices

    1172: #ifdef MOAB_HAVE_MPI
    1173:     /* filter all the non-owned and shared entities out of the list */
    1174:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
    1175:     MBERRNM(merr);
    1176:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT);
    1177:     MBERRNM(merr);
    1178: #endif

    1180:     /* get all the nodes via connectivity and the parent elements via adjacency information */
    1181:     merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false);
    1182:     MBERRNM(merr);
    1183:     merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION);
    1184:     MBERRNM(merr);
    1185:   } else {
    1186:     /* Let us query the hierarchy manager and get the results directly for this level */
    1187:     for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
    1188:       moab::EntityHandle elemHandle = *iter;
    1189:       if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
    1190:         dmmoab->bndyelems->insert(elemHandle);
    1191:         /* For this boundary element, query the vertices and add them to the list */
    1192:         std::vector<moab::EntityHandle> connect;
    1193:         merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect);
    1194:         MBERRNM(merr);
    1195:         for (unsigned iv = 0; iv < connect.size(); ++iv)
    1196:           if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv])) dmmoab->bndyvtx->insert(connect[iv]);
    1197:         /* Next, let us query the boundary faces and add them also to the list */
    1198:         std::vector<moab::EntityHandle> faces;
    1199:         merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim - 1, faces);
    1200:         MBERRNM(merr);
    1201:         for (unsigned ifa = 0; ifa < faces.size(); ++ifa)
    1202:           if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa])) dmmoab->bndyfaces->insert(faces[ifa]);
    1203:       }
    1204:     }
    1205: #ifdef MOAB_HAVE_MPI
    1206:     /* filter all the non-owned and shared entities out of the list */
    1207:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT);
    1208:     MBERRNM(merr);
    1209:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
    1210:     MBERRNM(merr);
    1211:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
    1212:     MBERRNM(merr);
    1213: #endif
    1214:   }
    1215:   PetscCall(PetscInfo(NULL, "Found %zu boundary vertices, %zu boundary faces and %zu boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size()));

    1217:   /* Get the material sets and populate the data for all locally owned elements */
    1218:   {
    1219:     PetscCall(PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials));
    1220:     /* Get the count of entities of particular type from dmmoab->elocal
    1221:        -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
    1222:     moab::Range msets;
    1223:     merr = dmmoab->mbiface->get_entities_by_type_and_tag(dmmoab->fileset, moab::MBENTITYSET, &dmmoab->material_tag, NULL, 1, msets, moab::Interface::UNION);
    1224:     MBERRNM(merr);
    1225:     if (msets.size() == 0) PetscCall(PetscInfo(NULL, "No material sets found in the fileset.\n"));

    1227:     for (unsigned i = 0; i < msets.size(); ++i) {
    1228:       moab::Range msetelems;
    1229:       merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);
    1230:       MBERRNM(merr);
    1231: #ifdef MOAB_HAVE_MPI
    1232:       /* filter all the non-owned and shared entities out of the list */
    1233:       merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
    1234:       MBERRNM(merr);
    1235: #endif

    1237:       int                partID;
    1238:       moab::EntityHandle mset = msets[i];
    1239:       merr                    = dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID);
    1240:       MBERRNM(merr);

    1242:       for (unsigned j = 0; j < msetelems.size(); ++j) dmmoab->materials[dmmoab->elocal->index(msetelems[j])] = partID;
    1243:     }
    1244:   }
    1245:   PetscFunctionReturn(PETSC_SUCCESS);
    1246: }

    1248: /*@C
    1249:   DMMoabCreateVertices - Creates and adds several vertices to the primary set represented by the DM.

    1251:   Collective

    1253:   Input Parameters:
    1254: + dm     - The `DM` object
    1255: . coords - The connectivity of the element
    1256: - nverts - The number of vertices that form the element

    1258:   Output Parameter:
    1259: . overts - The list of vertices that were created (can be `NULL`)

    1261:   Level: beginner

    1263: .seealso: `DMMOAB`, `DMMoabCreateSubmesh()`, `DMMoabCreateElement()`
    1264: @*/
    1265: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal *coords, PetscInt nverts, moab::Range *overts)
    1266: {
    1267:   moab::ErrorCode merr;
    1268:   DM_Moab        *dmmoab;
    1269:   moab::Range     verts;

    1271:   PetscFunctionBegin;
    1273:   PetscAssertPointer(coords, 2);

    1275:   dmmoab = (DM_Moab *)dm->data;

    1277:   /* Insert new points */
    1278:   merr = dmmoab->mbiface->create_vertices(&coords[0], nverts, verts);
    1279:   MBERRNM(merr);
    1280:   merr = dmmoab->mbiface->add_entities(dmmoab->fileset, verts);
    1281:   MBERRNM(merr);

    1283:   if (overts) *overts = verts;
    1284:   PetscFunctionReturn(PETSC_SUCCESS);
    1285: }

    1287: /*@C
    1288:   DMMoabCreateElement - Adds an element of specified type to the primary set represented by the DM.

    1290:   Collective

    1292:   Input Parameters:
    1293: + dm     - The DM object
    1294: . type   - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
    1295: . conn   - The connectivity of the element
    1296: - nverts - The number of vertices that form the element

    1298:   Output Parameter:
    1299: . oelem - The handle to the element created and added to the `DM` object

    1301:   Level: beginner

    1303: .seealso: `DMMOAB`, `DMMoabCreateSubmesh()`, `DMMoabCreateVertices()`
    1304: @*/
    1305: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle *conn, PetscInt nverts, moab::EntityHandle *oelem)
    1306: {
    1307:   moab::ErrorCode    merr;
    1308:   DM_Moab           *dmmoab;
    1309:   moab::EntityHandle elem;

    1311:   PetscFunctionBegin;
    1313:   PetscAssertPointer(conn, 3);

    1315:   dmmoab = (DM_Moab *)dm->data;

    1317:   /* Insert new element */
    1318:   merr = dmmoab->mbiface->create_element(type, conn, nverts, elem);
    1319:   MBERRNM(merr);
    1320:   merr = dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1);
    1321:   MBERRNM(merr);

    1323:   if (oelem) *oelem = elem;
    1324:   PetscFunctionReturn(PETSC_SUCCESS);
    1325: }

    1327: /*@C
    1328:   DMMoabCreateSubmesh - Creates a sub-`DM` object with a set that contains all vertices/elements of the parent
    1329:   in addition to providing support for dynamic mesh modifications. This is useful for AMR calculations to
    1330:   create a DM object on a refined level.

    1332:   Collective

    1334:   Input Parameters:
    1335: . dm - The `DM` object

    1337:   Output Parameter:
    1338: . newdm - The sub `DM` object with updated set information

    1340:   Level: advanced

    1342: .seealso: `DMMOAB`, `DMCreate()`, `DMMoabCreateVertices()`, `DMMoabCreateElement()`
    1343: @*/
    1344: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
    1345: {
    1346:   DM_Moab        *dmmoab;
    1347:   DM_Moab        *ndmmoab;
    1348:   moab::ErrorCode merr;

    1350:   PetscFunctionBegin;

    1353:   dmmoab = (DM_Moab *)dm->data;

    1355:   /* Create the basic DMMOAB object and keep the default parameters created by DM impls */
    1356:   PetscCall(DMMoabCreateMoab(((PetscObject)dm)->comm, dmmoab->mbiface, &dmmoab->ltog_tag, NULL, newdm));

    1358:   /* get all the necessary handles from the private DM object */
    1359:   ndmmoab = (DM_Moab *)(*newdm)->data;

    1361:   /* set the sub-mesh's parent DM reference */
    1362:   ndmmoab->parent = &dm;

    1364:   /* create a file set to associate all entities in current mesh */
    1365:   merr = ndmmoab->mbiface->create_meshset(moab::MESHSET_SET, ndmmoab->fileset);
    1366:   MBERR("Creating file set failed", merr);

    1368:   /* create a meshset and then add old fileset as child */
    1369:   merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal);
    1370:   MBERR("Adding child vertices to parent failed", merr);
    1371:   merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal);
    1372:   MBERR("Adding child elements to parent failed", merr);

    1374:   /* preserve the field association between the parent and sub-mesh objects */
    1375:   PetscCall(DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames));
    1376:   PetscFunctionReturn(PETSC_SUCCESS);
    1377: }

    1379: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
    1380: {
    1381:   DM_Moab    *dmmoab = (DM_Moab *)dm->data;
    1382:   const char *name;
    1383:   MPI_Comm    comm;
    1384:   PetscMPIInt size;

    1386:   PetscFunctionBegin;
    1387:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
    1388:   PetscCallMPI(MPI_Comm_size(comm, &size));
    1389:   PetscCall(PetscObjectGetName((PetscObject)dm, &name));
    1390:   PetscCall(PetscViewerASCIIPushTab(viewer));
    1391:   if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimensions:\n", name, dmmoab->dim));
    1392:   else PetscCall(PetscViewerASCIIPrintf(viewer, "Mesh in %" PetscInt_FMT " dimensions:\n", dmmoab->dim));
    1393:   /* print details about the global mesh */
    1394:   {
    1395:     PetscCall(PetscViewerASCIIPushTab(viewer));
    1396:     PetscCall(PetscViewerASCIIPrintf(viewer, "Sizes: cells=%" PetscInt_FMT ", vertices=%" PetscInt_FMT ", blocks=%" PetscInt_FMT "\n", dmmoab->nele, dmmoab->n, dmmoab->bs));
    1397:     /* print boundary data */
    1398:     PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary trace:\n"));
    1399:     {
    1400:       PetscCall(PetscViewerASCIIPushTab(viewer));
    1401:       PetscCall(PetscViewerASCIIPrintf(viewer, "cells=%zu, faces=%zu, vertices=%zu\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size()));
    1402:       PetscCall(PetscViewerASCIIPopTab(viewer));
    1403:     }
    1404:     /* print field data */
    1405:     PetscCall(PetscViewerASCIIPrintf(viewer, "Fields: %" PetscInt_FMT " components\n", dmmoab->numFields));
    1406:     {
    1407:       PetscCall(PetscViewerASCIIPushTab(viewer));
    1408:       for (int i = 0; i < dmmoab->numFields; ++i) PetscCall(PetscViewerASCIIPrintf(viewer, "[%" PetscInt_FMT "] - %s\n", i, dmmoab->fieldNames[i]));
    1409:       PetscCall(PetscViewerASCIIPopTab(viewer));
    1410:     }
    1411:     PetscCall(PetscViewerASCIIPopTab(viewer));
    1412:   }
    1413:   PetscCall(PetscViewerASCIIPopTab(viewer));
    1414:   PetscCall(PetscViewerFlush(viewer));
    1415:   PetscFunctionReturn(PETSC_SUCCESS);
    1416: }

    1418: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
    1419: {
    1420:   PetscFunctionReturn(PETSC_SUCCESS);
    1421: }

    1423: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
    1424: {
    1425:   PetscFunctionReturn(PETSC_SUCCESS);
    1426: }

    1428: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
    1429: {
    1430:   PetscBool iascii, ishdf5, isvtk;

    1432:   PetscFunctionBegin;
    1435:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
    1436:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
    1437:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
    1438:   if (iascii) {
    1439:     PetscCall(DMMoabView_Ascii(dm, viewer));
    1440:   } else if (ishdf5) {
    1441: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
    1442:     PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ));
    1443:     PetscCall(DMMoabView_HDF5(dm, viewer));
    1444:     PetscCall(PetscViewerPopFormat(viewer));
    1445: #else
    1446:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
    1447: #endif
    1448:   } else if (isvtk) {
    1449:     PetscCall(DMMoabView_VTK(dm, viewer));
    1450:   }
    1451:   PetscFunctionReturn(PETSC_SUCCESS);
    1452: }

    1454: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
    1455: {
    1456:   PetscFunctionBegin;
    1457:   dm->ops->view                     = DMView_Moab;
    1458:   dm->ops->load                     = NULL /* DMLoad_Moab */;
    1459:   dm->ops->setfromoptions           = DMSetFromOptions_Moab;
    1460:   dm->ops->clone                    = DMClone_Moab;
    1461:   dm->ops->setup                    = DMSetUp_Moab;
    1462:   dm->ops->createlocalsection       = NULL;
    1463:   dm->ops->createsectionpermutation = NULL;
    1464:   dm->ops->createdefaultconstraints = NULL;
    1465:   dm->ops->createglobalvector       = DMCreateGlobalVector_Moab;
    1466:   dm->ops->createlocalvector        = DMCreateLocalVector_Moab;
    1467:   dm->ops->getlocaltoglobalmapping  = NULL;
    1468:   dm->ops->createfieldis            = NULL;
    1469:   dm->ops->createcoordinatedm       = NULL /* DMCreateCoordinateDM_Moab */;
    1470:   dm->ops->getcoloring              = NULL;
    1471:   dm->ops->creatematrix             = DMCreateMatrix_Moab;
    1472:   dm->ops->createinterpolation      = DMCreateInterpolation_Moab;
    1473:   dm->ops->createinjection          = NULL /* DMCreateInjection_Moab */;
    1474:   dm->ops->refine                   = DMRefine_Moab;
    1475:   dm->ops->coarsen                  = DMCoarsen_Moab;
    1476:   dm->ops->refinehierarchy          = DMRefineHierarchy_Moab;
    1477:   dm->ops->coarsenhierarchy         = DMCoarsenHierarchy_Moab;
    1478:   dm->ops->globaltolocalbegin       = DMGlobalToLocalBegin_Moab;
    1479:   dm->ops->globaltolocalend         = DMGlobalToLocalEnd_Moab;
    1480:   dm->ops->localtoglobalbegin       = DMLocalToGlobalBegin_Moab;
    1481:   dm->ops->localtoglobalend         = DMLocalToGlobalEnd_Moab;
    1482:   dm->ops->destroy                  = DMDestroy_Moab;
    1483:   dm->ops->createsubdm              = NULL /* DMCreateSubDM_Moab */;
    1484:   dm->ops->getdimpoints             = NULL /* DMGetDimPoints_Moab */;
    1485:   dm->ops->locatepoints             = NULL /* DMLocatePoints_Moab */;
    1486:   PetscFunctionReturn(PETSC_SUCCESS);
    1487: }

    1489: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
    1490: {
    1491:   PetscFunctionBegin;
    1492:   /* get all the necessary handles from the private DM object */
    1493:   (*newdm)->data = (DM_Moab *)dm->data;
    1494:   ((DM_Moab *)dm->data)->refct++;

    1496:   PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMMOAB));
    1497:   PetscCall(DMInitialize_Moab(*newdm));
    1498:   PetscFunctionReturn(PETSC_SUCCESS);
    1499: }

    1501: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
    1502: {
    1503:   PetscFunctionBegin;
    1505:   PetscCall(PetscNew((DM_Moab **)&dm->data));

    1507:   ((DM_Moab *)dm->data)->bs            = 1;
    1508:   ((DM_Moab *)dm->data)->numFields     = 1;
    1509:   ((DM_Moab *)dm->data)->n             = 0;
    1510:   ((DM_Moab *)dm->data)->nloc          = 0;
    1511:   ((DM_Moab *)dm->data)->nghost        = 0;
    1512:   ((DM_Moab *)dm->data)->nele          = 0;
    1513:   ((DM_Moab *)dm->data)->neleloc       = 0;
    1514:   ((DM_Moab *)dm->data)->neleghost     = 0;
    1515:   ((DM_Moab *)dm->data)->ltog_map      = NULL;
    1516:   ((DM_Moab *)dm->data)->ltog_sendrecv = NULL;

    1518:   ((DM_Moab *)dm->data)->refct  = 1;
    1519:   ((DM_Moab *)dm->data)->parent = NULL;
    1520:   ((DM_Moab *)dm->data)->vlocal = new moab::Range();
    1521:   ((DM_Moab *)dm->data)->vowned = new moab::Range();
    1522:   ((DM_Moab *)dm->data)->vghost = new moab::Range();
    1523:   ((DM_Moab *)dm->data)->elocal = new moab::Range();
    1524:   ((DM_Moab *)dm->data)->eghost = new moab::Range();

    1526:   PetscCall(DMInitialize_Moab(dm));
    1527:   PetscFunctionReturn(PETSC_SUCCESS);
    1528: }