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.
 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:   Reference: https://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html

 19:   Level: intermediate

 21: .seealso: `DMType`, `DMMoabCreate()`, `DMCreate()`, `DMSetType()`, `DMMoabCreateMoab()`
 22: M*/

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

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

 52: /*@C
 53:   DMMoabCreate - Creates a DMMoab object, which encapsulates a moab instance

 55:   Collective

 57:   Input Parameter:
 58: . comm - The communicator for the DMMoab object

 60:   Output Parameter:
 61: . dmb  - The DMMoab object

 63:   Level: beginner

 65: @*/
 66: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
 67: {
 68:   PetscFunctionBegin;
 70:   PetscCall(DMCreate(comm, dmb));
 71:   PetscCall(DMSetType(*dmb, DMMOAB));
 72:   PetscFunctionReturn(PETSC_SUCCESS);
 73: }

 75: /*@C
 76:   DMMoabCreateMoab - Creates a DMMoab object, optionally from an instance and other data

 78:   Collective

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

 87:   Output Parameter:
 88: . dmb  - The DMMoab object

 90:   Level: intermediate

 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;

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: @*/
176: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
177: {
178:   PetscFunctionBegin;
180:   *pcomm = ((DM_Moab *)(dm)->data)->pcomm;
181:   PetscFunctionReturn(PETSC_SUCCESS);
182: }

184: #endif /* MOAB_HAVE_MPI */

186: /*@C
187:   DMMoabSetInterface - Set the MOAB instance used with this DMMoab

189:   Collective

191:   Input Parameters:
192: + dm      - The DMMoab object being set
193: - mbiface - The MOAB instance being set on this DMMoab

195:   Level: beginner

197: @*/
198: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
199: {
200:   DM_Moab *dmmoab = (DM_Moab *)(dm)->data;

202:   PetscFunctionBegin;
205: #ifdef MOAB_HAVE_MPI
206:   dmmoab->pcomm = NULL;
207: #endif
208:   dmmoab->mbiface          = mbiface;
209:   dmmoab->icreatedinstance = PETSC_FALSE;
210:   PetscFunctionReturn(PETSC_SUCCESS);
211: }

213: /*@C
214:   DMMoabGetInterface - Get the MOAB instance used with this DMMoab

216:   Collective

218:   Input Parameter:
219: . dm      - The DMMoab object being set

221:   Output Parameter:
222: . mbiface - The MOAB instance set on this DMMoab

224:   Level: beginner

226: @*/
227: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
228: {
229:   static PetscBool cite = PETSC_FALSE;

231:   PetscFunctionBegin;
233:   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, "
234:                                    "K. and Stimpson, C. and Ernst, C.},\n  year = {2004},  note = {Report}\n}\n",
235:                                    &cite));
236:   *mbiface = ((DM_Moab *)dm->data)->mbiface;
237:   PetscFunctionReturn(PETSC_SUCCESS);
238: }

240: /*@C
241:   DMMoabSetLocalVertices - Set the entities having DOFs on this DMMoab

243:   Collective

245:   Input Parameters:
246: + dm    - The DMMoab object being set
247: - range - The entities treated by this DMMoab

249:   Level: beginner

251: @*/
252: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
253: {
254:   moab::Range tmpvtxs;
255:   DM_Moab    *dmmoab = (DM_Moab *)(dm)->data;

257:   PetscFunctionBegin;
259:   dmmoab->vlocal->clear();
260:   dmmoab->vowned->clear();

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

264: #ifdef MOAB_HAVE_MPI
265:   moab::ErrorCode merr;
266:   /* filter based on parallel status */
267:   merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
268:   MBERRNM(merr);

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

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

291: /*@C
292:   DMMoabGetAllVertices - Get the entities having DOFs on this DMMoab

294:   Collective

296:   Input Parameter:
297: . dm    - The DMMoab object being set

299:   Output Parameter:
300: . owned - The local vertex entities in this DMMoab = (owned+ghosted)

302:   Level: beginner

304: @*/
305: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
306: {
307:   PetscFunctionBegin;
309:   if (local) *local = *((DM_Moab *)dm->data)->vlocal;
310:   PetscFunctionReturn(PETSC_SUCCESS);
311: }

313: /*@C
314:   DMMoabGetLocalVertices - Get the entities having DOFs on this DMMoab

316:   Collective

318:   Input Parameter:
319: . dm    - The DMMoab object being set

321:   Output Parameters:
322: + owned - The owned vertex entities in this DMMoab
323: - ghost - The ghosted entities (non-owned) stored locally in this partition

325:   Level: beginner

327: @*/
328: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
329: {
330:   PetscFunctionBegin;
332:   if (owned) *owned = ((DM_Moab *)dm->data)->vowned;
333:   if (ghost) *ghost = ((DM_Moab *)dm->data)->vghost;
334:   PetscFunctionReturn(PETSC_SUCCESS);
335: }

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

340:   Collective

342:   Input Parameter:
343: . dm    - The DMMoab object being set

345:   Output Parameter:
346: . range - The entities owned locally

348:   Level: beginner

350: @*/
351: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
352: {
353:   PetscFunctionBegin;
355:   if (range) *range = ((DM_Moab *)dm->data)->elocal;
356:   PetscFunctionReturn(PETSC_SUCCESS);
357: }

359: /*@C
360:   DMMoabSetLocalElements - Set the entities having DOFs on this DMMoab

362:   Collective

364:   Input Parameters:
365: + dm    - The DMMoab object being set
366: - range - The entities treated by this DMMoab

368:   Level: beginner

370: @*/
371: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
372: {
373:   DM_Moab *dmmoab = (DM_Moab *)(dm)->data;

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

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

400:   Collective

402:   Input Parameters:
403: + dm      - The DMMoab object being set
404: - ltogtag - The MOAB tag used for local to global ids

406:   Level: beginner

408: @*/
409: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
410: {
411:   PetscFunctionBegin;
413:   ((DM_Moab *)dm->data)->ltog_tag = ltogtag;
414:   PetscFunctionReturn(PETSC_SUCCESS);
415: }

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

420:   Collective

422:   Input Parameter:
423: . dm      - The DMMoab object being set

425:   Output Parameter:
426: . ltogtag - The MOAB tag used for local to global ids

428:   Level: beginner

430: @*/
431: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
432: {
433:   PetscFunctionBegin;
435:   *ltog_tag = ((DM_Moab *)dm->data)->ltog_tag;
436:   PetscFunctionReturn(PETSC_SUCCESS);
437: }

439: /*@C
440:   DMMoabSetBlockSize - Set the block size used with this DMMoab

442:   Collective

444:   Input Parameters:
445: + dm - The DMMoab object being set
446: - bs - The block size used with this DMMoab

448:   Level: beginner

450: @*/
451: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
452: {
453:   PetscFunctionBegin;
455:   ((DM_Moab *)dm->data)->bs = bs;
456:   PetscFunctionReturn(PETSC_SUCCESS);
457: }

459: /*@C
460:   DMMoabGetBlockSize - Get the block size used with this DMMoab

462:   Collective

464:   Input Parameter:
465: . dm - The DMMoab object being set

467:   Output Parameter:
468: . bs - The block size used with this DMMoab

470:   Level: beginner

472: @*/
473: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
474: {
475:   PetscFunctionBegin;
477:   *bs = ((DM_Moab *)dm->data)->bs;
478:   PetscFunctionReturn(PETSC_SUCCESS);
479: }

481: /*@C
482:   DMMoabGetSize - Get the global vertex size used with this DMMoab

484:   Collective on dm

486:   Input Parameter:
487: . dm - The DMMoab object being set

489:   Output Parameters:
490: + neg - The number of global elements in the DMMoab instance
491: - nvg - The number of global vertices in the DMMoab instance

493:   Level: beginner

495: @*/
496: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
497: {
498:   PetscFunctionBegin;
500:   if (neg) *neg = ((DM_Moab *)dm->data)->nele;
501:   if (nvg) *nvg = ((DM_Moab *)dm->data)->n;
502:   PetscFunctionReturn(PETSC_SUCCESS);
503: }

505: /*@C
506:   DMMoabGetLocalSize - Get the local and ghosted vertex size used with this DMMoab

508:   Collective on dm

510:   Input Parameter:
511: . dm - The DMMoab object being set

513:   Output Parameters:
514: + nel - The number of owned elements in this processor
515: . neg - The number of ghosted elements in this processor
516: . nvl - The number of owned vertices in this processor
517: - nvg - The number of ghosted vertices in this processor

519:   Level: beginner

521: @*/
522: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
523: {
524:   PetscFunctionBegin;
526:   if (nel) *nel = ((DM_Moab *)dm->data)->neleloc;
527:   if (neg) *neg = ((DM_Moab *)dm->data)->neleghost;
528:   if (nvl) *nvl = ((DM_Moab *)dm->data)->nloc;
529:   if (nvg) *nvg = ((DM_Moab *)dm->data)->nghost;
530:   PetscFunctionReturn(PETSC_SUCCESS);
531: }

533: /*@C
534:   DMMoabGetOffset - Get the local offset for the global vector

536:   Collective

538:   Input Parameter:
539: . dm - The DMMoab object being set

541:   Output Parameter:
542: . offset - The local offset for the global vector

544:   Level: beginner

546: @*/
547: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
548: {
549:   PetscFunctionBegin;
551:   *offset = ((DM_Moab *)dm->data)->vstart;
552:   PetscFunctionReturn(PETSC_SUCCESS);
553: }

555: /*@C
556:   DMMoabGetDimension - Get the dimension of the DM Mesh

558:   Collective

560:   Input Parameter:
561: . dm - The DMMoab object

563:   Output Parameter:
564: . dim - The dimension of DM

566:   Level: beginner

568: @*/
569: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
570: {
571:   PetscFunctionBegin;
573:   *dim = ((DM_Moab *)dm->data)->dim;
574:   PetscFunctionReturn(PETSC_SUCCESS);
575: }

577: /*@C
578:   DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
579:   generated through uniform refinement.

581:   Collective on dm

583:   Input Parameter:
584: . dm - The DMMoab object being set

586:   Output Parameter:
587: . nvg - The current mesh hierarchy level

589:   Level: beginner

591: @*/
592: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
593: {
594:   PetscFunctionBegin;
596:   if (nlevel) *nlevel = ((DM_Moab *)dm->data)->hlevel;
597:   PetscFunctionReturn(PETSC_SUCCESS);
598: }

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

603:   Collective

605:   Input Parameters:
606: + dm - The DMMoab object
607: - ehandle - The element entity handle

609:   Output Parameter:
610: . mat - The material ID for the current entity

612:   Level: beginner

614: @*/
615: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
616: {
617:   DM_Moab *dmmoab;

619:   PetscFunctionBegin;
621:   if (*mat) {
622:     dmmoab = (DM_Moab *)(dm)->data;
623:     *mat   = dmmoab->materials[dmmoab->elocal->index(ehandle)];
624:   }
625:   PetscFunctionReturn(PETSC_SUCCESS);
626: }

628: /*@C
629:   DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities

631:   Collective

633:   Input Parameters:
634: + dm - The DMMoab object
635: . nconn - Number of entities whose coordinates are needed
636: - conn - The vertex entity handles

638:   Output Parameter:
639: . vpos - The coordinates of the requested vertex entities

641:   Level: beginner

643: .seealso: `DMMoabGetVertexConnectivity()`
644: @*/
645: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
646: {
647:   DM_Moab        *dmmoab;
648:   moab::ErrorCode merr;

650:   PetscFunctionBegin;
654:   dmmoab = (DM_Moab *)(dm)->data;

656:   /* Get connectivity information in MOAB canonical ordering */
657:   if (dmmoab->hlevel) {
658:     merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle *>(conn), nconn, dmmoab->hlevel, vpos);
659:     MBERRNM(merr);
660:   } else {
661:     merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);
662:     MBERRNM(merr);
663:   }
664:   PetscFunctionReturn(PETSC_SUCCESS);
665: }

667: /*@C
668:   DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity

670:   Collective

672:   Input Parameters:
673: + dm - The DMMoab object
674: - vhandle - Vertex entity handle

676:   Output Parameters:
677: + nconn - Number of entities whose coordinates are needed
678: - conn - The vertex entity handles

680:   Level: beginner

682: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabRestoreVertexConnectivity()`
683: @*/
684: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt *nconn, moab::EntityHandle **conn)
685: {
686:   DM_Moab                        *dmmoab;
687:   std::vector<moab::EntityHandle> adj_entities, connect;
688:   moab::ErrorCode                 merr;

690:   PetscFunctionBegin;
693:   dmmoab = (DM_Moab *)(dm)->data;

695:   /* Get connectivity information in MOAB canonical ordering */
696:   merr = dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION);
697:   MBERRNM(merr);
698:   merr = dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect);
699:   MBERRNM(merr);

701:   if (conn) {
702:     PetscCall(PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn));
703:     PetscCall(PetscArraycpy(*conn, &connect[0], connect.size()));
704:   }
705:   if (nconn) *nconn = connect.size();
706:   PetscFunctionReturn(PETSC_SUCCESS);
707: }

709: /*@C
710:   DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity

712:   Collective

714:   Input Parameters:
715: + dm - The DMMoab object
716: . vhandle - Vertex entity handle
717: . nconn - Number of entities whose coordinates are needed
718: - conn - The vertex entity handles

720:   Level: beginner

722: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`
723: @*/
724: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, moab::EntityHandle **conn)
725: {
726:   PetscFunctionBegin;

730:   if (conn) PetscCall(PetscFree(*conn));
731:   if (nconn) *nconn = 0;
732:   PetscFunctionReturn(PETSC_SUCCESS);
733: }

735: /*@C
736:   DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity

738:   Collective

740:   Input Parameters:
741: + dm - The DMMoab object
742: - ehandle - Vertex entity handle

744:   Output Parameters:
745: + nconn - Number of entities whose coordinates are needed
746: - conn - The vertex entity handles

748:   Level: beginner

750: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`, `DMMoabRestoreVertexConnectivity()`
751: @*/
752: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, const moab::EntityHandle **conn)
753: {
754:   DM_Moab                        *dmmoab;
755:   const moab::EntityHandle       *connect;
756:   std::vector<moab::EntityHandle> vconn;
757:   moab::ErrorCode                 merr;
758:   PetscInt                        nnodes;

760:   PetscFunctionBegin;
763:   dmmoab = (DM_Moab *)(dm)->data;

765:   /* Get connectivity information in MOAB canonical ordering */
766:   merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes);
767:   MBERRNM(merr);
768:   if (conn) *conn = connect;
769:   if (nconn) *nconn = nnodes;
770:   PetscFunctionReturn(PETSC_SUCCESS);
771: }

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

776:   Collective

778:   Input Parameters:
779: + dm - The DMMoab object
780: - ent - Entity handle

782:   Output Parameter:
783: . ent_on_boundary - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise

785:   Level: beginner

787: .seealso: `DMMoabCheckBoundaryVertices()`
788: @*/
789: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool *ent_on_boundary)
790: {
791:   moab::EntityType etype;
792:   DM_Moab         *dmmoab;
793:   PetscInt         edim;

795:   PetscFunctionBegin;
798:   dmmoab = (DM_Moab *)(dm)->data;

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

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

807:   *ent_on_boundary = PETSC_FALSE;
808:   if (etype == moab::MBVERTEX && edim == 0) {
809:     *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
810:   } else {
811:     if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
812:       if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
813:     } else { /* next check the lower-dimensional faces */
814:       if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
815:     }
816:   }
817:   PetscFunctionReturn(PETSC_SUCCESS);
818: }

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

823:   Input Parameters:
824: + dm - The DMMoab object
825: . nconn - Number of handles
826: - cnt - Array of entity handles

828:   Output Parameter:
829: . isbdvtx - Array of boundary markers - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise

831:   Level: beginner

833: .seealso: `DMMoabIsEntityOnBoundary()`
834: @*/
835: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool *isbdvtx)
836: {
837:   DM_Moab *dmmoab;
838:   PetscInt i;

840:   PetscFunctionBegin;
844:   dmmoab = (DM_Moab *)(dm)->data;

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

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

853:   Input Parameter:
854: . dm - The DMMoab object

856:   Output Parameters:
857: + bdvtx - Boundary vertices
858: . bdelems - Boundary elements
859: - bdfaces - Boundary faces

861:   Level: beginner

863: .seealso: `DMMoabCheckBoundaryVertices()`, `DMMoabIsEntityOnBoundary()`
864: @*/
865: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range **bdelems, const moab::Range **bdfaces)
866: {
867:   DM_Moab *dmmoab;

869:   PetscFunctionBegin;
871:   dmmoab = (DM_Moab *)(dm)->data;

873:   if (bdvtx) *bdvtx = dmmoab->bndyvtx;
874:   if (bdfaces) *bdfaces = dmmoab->bndyfaces;
875:   if (bdelems) *bdfaces = dmmoab->bndyelems;
876:   PetscFunctionReturn(PETSC_SUCCESS);
877: }

879: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
880: {
881:   PetscInt        i;
882:   moab::ErrorCode merr;
883:   DM_Moab        *dmmoab = (DM_Moab *)dm->data;

885:   PetscFunctionBegin;

888:   dmmoab->refct--;
889:   if (!dmmoab->refct) {
890:     delete dmmoab->vlocal;
891:     delete dmmoab->vowned;
892:     delete dmmoab->vghost;
893:     delete dmmoab->elocal;
894:     delete dmmoab->eghost;
895:     delete dmmoab->bndyvtx;
896:     delete dmmoab->bndyfaces;
897:     delete dmmoab->bndyelems;

899:     PetscCall(PetscFree(dmmoab->gsindices));
900:     PetscCall(PetscFree2(dmmoab->gidmap, dmmoab->lidmap));
901:     PetscCall(PetscFree(dmmoab->dfill));
902:     PetscCall(PetscFree(dmmoab->ofill));
903:     PetscCall(PetscFree(dmmoab->materials));
904:     if (dmmoab->fieldNames) {
905:       for (i = 0; i < dmmoab->numFields; i++) PetscCall(PetscFree(dmmoab->fieldNames[i]));
906:       PetscCall(PetscFree(dmmoab->fieldNames));
907:     }

909:     if (dmmoab->nhlevels) {
910:       PetscCall(PetscFree(dmmoab->hsets));
911:       dmmoab->nhlevels = 0;
912:       if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
913:       dmmoab->hierarchy = NULL;
914:     }

916:     if (dmmoab->icreatedinstance) {
917:       delete dmmoab->pcomm;
918:       merr = dmmoab->mbiface->delete_mesh();
919:       MBERRNM(merr);
920:       delete dmmoab->mbiface;
921:     }
922:     dmmoab->mbiface = NULL;
923: #ifdef MOAB_HAVE_MPI
924:     dmmoab->pcomm = NULL;
925: #endif
926:     PetscCall(VecScatterDestroy(&dmmoab->ltog_sendrecv));
927:     PetscCall(ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map));
928:     PetscCall(PetscFree(dm->data));
929:   }
930:   PetscFunctionReturn(PETSC_SUCCESS);
931: }

933: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(DM dm, PetscOptionItems *PetscOptionsObject)
934: {
935:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

937:   PetscFunctionBegin;
938:   PetscOptionsHeadBegin(PetscOptionsObject, "DMMoab Options");
939:   PetscCall(PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL, 0));
940:   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));
941:   /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
942:   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));
943:   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));
944:   PetscCall(PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum *)&dmmoab->read_mode, NULL));
945:   PetscCall(PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum *)&dmmoab->write_mode, NULL));
946:   PetscOptionsHeadEnd();
947:   PetscFunctionReturn(PETSC_SUCCESS);
948: }

950: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
951: {
952:   moab::ErrorCode       merr;
953:   Vec                   local, global;
954:   IS                    from, to;
955:   moab::Range::iterator iter;
956:   PetscInt              i, j, f, bs, vent, totsize, *lgmap;
957:   DM_Moab              *dmmoab = (DM_Moab *)dm->data;
958:   moab::Range           adjs;

960:   PetscFunctionBegin;
962:   /* Get the local and shared vertices and cache it */
963:   PetscCheck(dmmoab->mbiface != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB Interface before calling SetUp.");
964: #ifdef MOAB_HAVE_MPI
965:   PetscCheck(dmmoab->pcomm != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB ParallelComm object before calling SetUp.");
966: #endif

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

974: #ifdef MOAB_HAVE_MPI
975:     /* filter based on parallel status */
976:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
977:     MBERRNM(merr);

979:     /* filter all the non-owned and shared entities out of the list */
980:     // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
981:     adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
982:     merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
983:     MBERRNM(merr);
984:     adjs            = moab::subtract(adjs, *dmmoab->vghost);
985:     *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
986: #else
987:     *dmmoab->vowned = *dmmoab->vlocal;
988: #endif

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

994: #ifdef MOAB_HAVE_MPI
995:     PetscCall(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
996:     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));
997: #else
998:     dmmoab->n       = dmmoab->nloc;
999: #endif
1000:   }

1002:   {
1003:     /* get the information about the local elements in the mesh */
1004:     dmmoab->eghost->clear();

1006:     /* first decipher the leading dimension */
1007:     for (i = 3; i > 0; i--) {
1008:       dmmoab->elocal->clear();
1009:       merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false);
1010:       MBERRNM(merr);

1012:       /* store the current mesh dimension */
1013:       if (dmmoab->elocal->size()) {
1014:         dmmoab->dim = i;
1015:         break;
1016:       }
1017:     }

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

1021: #ifdef MOAB_HAVE_MPI
1022:     /* filter the ghosted and owned element list */
1023:     *dmmoab->eghost = *dmmoab->elocal;
1024:     merr            = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1025:     MBERRNM(merr);
1026:     *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
1027: #endif

1029:     dmmoab->neleloc   = dmmoab->elocal->size();
1030:     dmmoab->neleghost = dmmoab->eghost->size();

1032: #ifdef MOAB_HAVE_MPI
1033:     PetscCall(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
1034:     PetscCall(PetscInfo(NULL, "%d-dim elements: owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost));
1035: #else
1036:     dmmoab->nele    = dmmoab->neleloc;
1037: #endif
1038:   }

1040:   bs = dmmoab->bs;
1041:   if (!dmmoab->ltog_tag) {
1042:     /* Get the global ID tag. The global ID tag is applied to each
1043:        vertex. It acts as an global identifier which MOAB uses to
1044:        assemble the individual pieces of the mesh */
1045:     merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
1046:     MBERRNM(merr);
1047:   }

1049:   totsize = dmmoab->vlocal->size();
1050:   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);
1051:   PetscCall(PetscCalloc1(totsize, &dmmoab->gsindices));
1052:   {
1053:     /* first get the local indices */
1054:     merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]);
1055:     MBERRNM(merr);
1056:     if (dmmoab->nghost) { /* next get the ghosted indices */
1057:       merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]);
1058:       MBERRNM(merr);
1059:     }

1061:     /* find out the local and global minima of GLOBAL_ID */
1062:     dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1063:     for (i = 0; i < totsize; ++i) {
1064:       if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1065:       if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1066:     }

1068:     PetscCall(MPIU_Allreduce(&dmmoab->lminmax[0], &dmmoab->gminmax[0], 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm));
1069:     PetscCall(MPIU_Allreduce(&dmmoab->lminmax[1], &dmmoab->gminmax[1], 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm));

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

1076:     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]));
1077:   }
1078:   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,
1079:              dmmoab->numFields);

1081:   {
1082:     dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1083:     dmmoab->seqend   = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1084:     PetscCall(PetscInfo(NULL, "SEQUENCE: Local [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "]\n", dmmoab->seqstart, dmmoab->seqend));

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

1089:     i = j = 0;
1090:     /* set the owned vertex data first */
1091:     for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1092:       vent                 = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1093:       dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1094:       dmmoab->lidmap[vent] = i;
1095:       for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1096:     }
1097:     /* next arrange all the ghosted data information */
1098:     for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1099:       vent                 = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1100:       dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1101:       dmmoab->lidmap[vent] = i;
1102:       for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1103:     }

1105:     /* We need to create the Global to Local Vector Scatter Contexts
1106:        1) First create a local and global vector
1107:        2) Create a local and global IS
1108:        3) Create VecScatter and LtoGMapping objects
1109:        4) Cleanup the IS and Vec objects
1110:     */
1111:     PetscCall(DMCreateGlobalVector(dm, &global));
1112:     PetscCall(DMCreateLocalVector(dm, &local));

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

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

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

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

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

1131:     /* clean up IS, Vec */
1132:     PetscCall(PetscFree(lgmap));
1133:     PetscCall(ISDestroy(&from));
1134:     PetscCall(ISDestroy(&to));
1135:     PetscCall(VecDestroy(&local));
1136:     PetscCall(VecDestroy(&global));
1137:   }

1139:   dmmoab->bndyvtx   = new moab::Range();
1140:   dmmoab->bndyfaces = new moab::Range();
1141:   dmmoab->bndyelems = new moab::Range();
1142:   /* skin the boundary and store nodes */
1143:   if (!dmmoab->hlevel) {
1144:     /* get the skin vertices of boundary faces for the current partition and then filter
1145:        the local, boundary faces, vertices and elements alone via PSTATUS flags;
1146:        this should not give us any ghosted boundary, but if user needs such a functionality
1147:        it would be easy to add it based on the find_skin query below */
1148:     moab::Skinner skinner(dmmoab->mbiface);

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

1154: #ifdef MOAB_HAVE_MPI
1155:     /* filter all the non-owned and shared entities out of the list */
1156:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1157:     MBERRNM(merr);
1158:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT);
1159:     MBERRNM(merr);
1160: #endif

1162:     /* get all the nodes via connectivity and the parent elements via adjacency information */
1163:     merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false);
1164:     MBERRNM(merr);
1165:     merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION);
1166:     MBERRNM(merr);
1167:   } else {
1168:     /* Let us query the hierarchy manager and get the results directly for this level */
1169:     for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1170:       moab::EntityHandle elemHandle = *iter;
1171:       if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1172:         dmmoab->bndyelems->insert(elemHandle);
1173:         /* For this boundary element, query the vertices and add them to the list */
1174:         std::vector<moab::EntityHandle> connect;
1175:         merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect);
1176:         MBERRNM(merr);
1177:         for (unsigned iv = 0; iv < connect.size(); ++iv)
1178:           if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv])) dmmoab->bndyvtx->insert(connect[iv]);
1179:         /* Next, let us query the boundary faces and add them also to the list */
1180:         std::vector<moab::EntityHandle> faces;
1181:         merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim - 1, faces);
1182:         MBERRNM(merr);
1183:         for (unsigned ifa = 0; ifa < faces.size(); ++ifa)
1184:           if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa])) dmmoab->bndyfaces->insert(faces[ifa]);
1185:       }
1186:     }
1187: #ifdef MOAB_HAVE_MPI
1188:     /* filter all the non-owned and shared entities out of the list */
1189:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1190:     MBERRNM(merr);
1191:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1192:     MBERRNM(merr);
1193:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1194:     MBERRNM(merr);
1195: #endif
1196:   }
1197:   PetscCall(PetscInfo(NULL, "Found %zu boundary vertices, %zu boundary faces and %zu boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size()));

1199:   /* Get the material sets and populate the data for all locally owned elements */
1200:   {
1201:     PetscCall(PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials));
1202:     /* Get the count of entities of particular type from dmmoab->elocal
1203:        -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1204:     moab::Range msets;
1205:     merr = dmmoab->mbiface->get_entities_by_type_and_tag(dmmoab->fileset, moab::MBENTITYSET, &dmmoab->material_tag, NULL, 1, msets, moab::Interface::UNION);
1206:     MBERRNM(merr);
1207:     if (msets.size() == 0) PetscCall(PetscInfo(NULL, "No material sets found in the fileset."));

1209:     for (unsigned i = 0; i < msets.size(); ++i) {
1210:       moab::Range msetelems;
1211:       merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);
1212:       MBERRNM(merr);
1213: #ifdef MOAB_HAVE_MPI
1214:       /* filter all the non-owned and shared entities out of the list */
1215:       merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1216:       MBERRNM(merr);
1217: #endif

1219:       int                partID;
1220:       moab::EntityHandle mset = msets[i];
1221:       merr                    = dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID);
1222:       MBERRNM(merr);

1224:       for (unsigned j = 0; j < msetelems.size(); ++j) dmmoab->materials[dmmoab->elocal->index(msetelems[j])] = partID;
1225:     }
1226:   }

1228:   PetscFunctionReturn(PETSC_SUCCESS);
1229: }

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

1234:   Collective

1236:   Input Parameters:
1237: + dm - The DM object
1238: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1239: . conn - The connectivity of the element
1240: - nverts - The number of vertices that form the element

1242:   Output Parameter:
1243: . overts  - The list of vertices that were created (can be NULL)

1245:   Level: beginner

1247: .seealso: `DMMoabCreateSubmesh()`, `DMMoabCreateElement()`
1248: @*/
1249: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal *coords, PetscInt nverts, moab::Range *overts)
1250: {
1251:   moab::ErrorCode merr;
1252:   DM_Moab        *dmmoab;
1253:   moab::Range     verts;

1255:   PetscFunctionBegin;

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

1261:   /* Insert new points */
1262:   merr = dmmoab->mbiface->create_vertices(&coords[0], nverts, verts);
1263:   MBERRNM(merr);
1264:   merr = dmmoab->mbiface->add_entities(dmmoab->fileset, verts);
1265:   MBERRNM(merr);

1267:   if (overts) *overts = verts;
1268:   PetscFunctionReturn(PETSC_SUCCESS);
1269: }

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

1274:   Collective

1276:   Input Parameters:
1277: + dm - The DM object
1278: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1279: . conn - The connectivity of the element
1280: - nverts - The number of vertices that form the element

1282:   Output Parameter:
1283: . oelem  - The handle to the element created and added to the DM object

1285:   Level: beginner

1287: .seealso: `DMMoabCreateSubmesh()`, `DMMoabCreateVertices()`
1288: @*/
1289: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle *conn, PetscInt nverts, moab::EntityHandle *oelem)
1290: {
1291:   moab::ErrorCode    merr;
1292:   DM_Moab           *dmmoab;
1293:   moab::EntityHandle elem;

1295:   PetscFunctionBegin;

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

1301:   /* Insert new element */
1302:   merr = dmmoab->mbiface->create_element(type, conn, nverts, elem);
1303:   MBERRNM(merr);
1304:   merr = dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1);
1305:   MBERRNM(merr);

1307:   if (oelem) *oelem = elem;
1308:   PetscFunctionReturn(PETSC_SUCCESS);
1309: }

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

1316:   Collective

1318:   Input Parameters:
1319: . dm - The DM object

1321:   Output Parameter:
1322: . newdm  - The sub DM object with updated set information

1324:   Level: advanced

1326: .seealso: `DMCreate()`, `DMMoabCreateVertices()`, `DMMoabCreateElement()`
1327: @*/
1328: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1329: {
1330:   DM_Moab        *dmmoab;
1331:   DM_Moab        *ndmmoab;
1332:   moab::ErrorCode merr;

1334:   PetscFunctionBegin;

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

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

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

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

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

1352:   /* create a meshset and then add old fileset as child */
1353:   merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal);
1354:   MBERR("Adding child vertices to parent failed", merr);
1355:   merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal);
1356:   MBERR("Adding child elements to parent failed", merr);

1358:   /* preserve the field association between the parent and sub-mesh objects */
1359:   PetscCall(DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames));
1360:   PetscFunctionReturn(PETSC_SUCCESS);
1361: }

1363: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1364: {
1365:   DM_Moab    *dmmoab = (DM_Moab *)(dm)->data;
1366:   const char *name;
1367:   MPI_Comm    comm;
1368:   PetscMPIInt size;

1370:   PetscFunctionBegin;
1371:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1372:   PetscCallMPI(MPI_Comm_size(comm, &size));
1373:   PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1374:   PetscCall(PetscViewerASCIIPushTab(viewer));
1375:   if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimensions:\n", name, dmmoab->dim));
1376:   else PetscCall(PetscViewerASCIIPrintf(viewer, "Mesh in %" PetscInt_FMT " dimensions:\n", dmmoab->dim));
1377:   /* print details about the global mesh */
1378:   {
1379:     PetscCall(PetscViewerASCIIPushTab(viewer));
1380:     PetscCall(PetscViewerASCIIPrintf(viewer, "Sizes: cells=%" PetscInt_FMT ", vertices=%" PetscInt_FMT ", blocks=%" PetscInt_FMT "\n", dmmoab->nele, dmmoab->n, dmmoab->bs));
1381:     /* print boundary data */
1382:     PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary trace:\n"));
1383:     {
1384:       PetscCall(PetscViewerASCIIPushTab(viewer));
1385:       PetscCall(PetscViewerASCIIPrintf(viewer, "cells=%zu, faces=%zu, vertices=%zu\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size()));
1386:       PetscCall(PetscViewerASCIIPopTab(viewer));
1387:     }
1388:     /* print field data */
1389:     PetscCall(PetscViewerASCIIPrintf(viewer, "Fields: %" PetscInt_FMT " components\n", dmmoab->numFields));
1390:     {
1391:       PetscCall(PetscViewerASCIIPushTab(viewer));
1392:       for (int i = 0; i < dmmoab->numFields; ++i) PetscCall(PetscViewerASCIIPrintf(viewer, "[%" PetscInt_FMT "] - %s\n", i, dmmoab->fieldNames[i]));
1393:       PetscCall(PetscViewerASCIIPopTab(viewer));
1394:     }
1395:     PetscCall(PetscViewerASCIIPopTab(viewer));
1396:   }
1397:   PetscCall(PetscViewerASCIIPopTab(viewer));
1398:   PetscCall(PetscViewerFlush(viewer));
1399:   PetscFunctionReturn(PETSC_SUCCESS);
1400: }

1402: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1403: {
1404:   PetscFunctionReturn(PETSC_SUCCESS);
1405: }

1407: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1408: {
1409:   PetscFunctionReturn(PETSC_SUCCESS);
1410: }

1412: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1413: {
1414:   PetscBool iascii, ishdf5, isvtk;

1416:   PetscFunctionBegin;
1419:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1420:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1421:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1422:   if (iascii) {
1423:     PetscCall(DMMoabView_Ascii(dm, viewer));
1424:   } else if (ishdf5) {
1425: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1426:     PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ));
1427:     PetscCall(DMMoabView_HDF5(dm, viewer));
1428:     PetscCall(PetscViewerPopFormat(viewer));
1429: #else
1430:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1431: #endif
1432:   } else if (isvtk) {
1433:     PetscCall(DMMoabView_VTK(dm, viewer));
1434:   }
1435:   PetscFunctionReturn(PETSC_SUCCESS);
1436: }

1438: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
1439: {
1440:   PetscFunctionBegin;
1441:   dm->ops->view                     = DMView_Moab;
1442:   dm->ops->load                     = NULL /* DMLoad_Moab */;
1443:   dm->ops->setfromoptions           = DMSetFromOptions_Moab;
1444:   dm->ops->clone                    = DMClone_Moab;
1445:   dm->ops->setup                    = DMSetUp_Moab;
1446:   dm->ops->createlocalsection       = NULL;
1447:   dm->ops->createdefaultconstraints = NULL;
1448:   dm->ops->createglobalvector       = DMCreateGlobalVector_Moab;
1449:   dm->ops->createlocalvector        = DMCreateLocalVector_Moab;
1450:   dm->ops->getlocaltoglobalmapping  = NULL;
1451:   dm->ops->createfieldis            = NULL;
1452:   dm->ops->createcoordinatedm       = NULL /* DMCreateCoordinateDM_Moab */;
1453:   dm->ops->getcoloring              = NULL;
1454:   dm->ops->creatematrix             = DMCreateMatrix_Moab;
1455:   dm->ops->createinterpolation      = DMCreateInterpolation_Moab;
1456:   dm->ops->createinjection          = NULL /* DMCreateInjection_Moab */;
1457:   dm->ops->refine                   = DMRefine_Moab;
1458:   dm->ops->coarsen                  = DMCoarsen_Moab;
1459:   dm->ops->refinehierarchy          = DMRefineHierarchy_Moab;
1460:   dm->ops->coarsenhierarchy         = DMCoarsenHierarchy_Moab;
1461:   dm->ops->globaltolocalbegin       = DMGlobalToLocalBegin_Moab;
1462:   dm->ops->globaltolocalend         = DMGlobalToLocalEnd_Moab;
1463:   dm->ops->localtoglobalbegin       = DMLocalToGlobalBegin_Moab;
1464:   dm->ops->localtoglobalend         = DMLocalToGlobalEnd_Moab;
1465:   dm->ops->destroy                  = DMDestroy_Moab;
1466:   dm->ops->createsubdm              = NULL /* DMCreateSubDM_Moab */;
1467:   dm->ops->getdimpoints             = NULL /* DMGetDimPoints_Moab */;
1468:   dm->ops->locatepoints             = NULL /* DMLocatePoints_Moab */;
1469:   PetscFunctionReturn(PETSC_SUCCESS);
1470: }

1472: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1473: {
1474:   PetscFunctionBegin;
1475:   /* get all the necessary handles from the private DM object */
1476:   (*newdm)->data = (DM_Moab *)dm->data;
1477:   ((DM_Moab *)dm->data)->refct++;

1479:   PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMMOAB));
1480:   PetscCall(DMInitialize_Moab(*newdm));
1481:   PetscFunctionReturn(PETSC_SUCCESS);
1482: }

1484: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1485: {
1486:   PetscFunctionBegin;
1488:   PetscCall(PetscNew((DM_Moab **)&dm->data));

1490:   ((DM_Moab *)dm->data)->bs            = 1;
1491:   ((DM_Moab *)dm->data)->numFields     = 1;
1492:   ((DM_Moab *)dm->data)->n             = 0;
1493:   ((DM_Moab *)dm->data)->nloc          = 0;
1494:   ((DM_Moab *)dm->data)->nghost        = 0;
1495:   ((DM_Moab *)dm->data)->nele          = 0;
1496:   ((DM_Moab *)dm->data)->neleloc       = 0;
1497:   ((DM_Moab *)dm->data)->neleghost     = 0;
1498:   ((DM_Moab *)dm->data)->ltog_map      = NULL;
1499:   ((DM_Moab *)dm->data)->ltog_sendrecv = NULL;

1501:   ((DM_Moab *)dm->data)->refct  = 1;
1502:   ((DM_Moab *)dm->data)->parent = NULL;
1503:   ((DM_Moab *)dm->data)->vlocal = new moab::Range();
1504:   ((DM_Moab *)dm->data)->vowned = new moab::Range();
1505:   ((DM_Moab *)dm->data)->vghost = new moab::Range();
1506:   ((DM_Moab *)dm->data)->elocal = new moab::Range();
1507:   ((DM_Moab *)dm->data)->eghost = new moab::Range();

1509:   PetscCall(DMInitialize_Moab(dm));
1510:   PetscFunctionReturn(PETSC_SUCCESS);
1511: }