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: {
 69:   DMCreate(comm, dmb);
 70:   DMSetType(*dmb, DMMOAB);
 71:   return 0;
 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 MOAB 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: @*/
 92: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
 93: {
 94:   moab::ErrorCode merr;
 95:   DM              dmmb;
 96:   DM_Moab        *dmmoab;


100:   DMMoabCreate(comm, &dmmb);
101:   dmmoab = (DM_Moab *)(dmmb)->data;

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

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

116: #ifdef MOAB_HAVE_MPI
117:   moab::EntityHandle partnset;

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

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

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

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

149:   merr = dmmoab->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dmmoab->material_tag);
150:   MBERRNM(merr);

152:   /* set the local range of entities (vertices) of interest */
153:   if (range) DMMoabSetLocalVertices(dmmb, range);
154:   *dmb = dmmb;
155:   return 0;
156: }

158: #ifdef MOAB_HAVE_MPI

160: /*@C
161:   DMMoabGetParallelComm - Get the ParallelComm used with this DMMoab

163:   Collective

165:   Input Parameter:
166: . dm    - The DMMoab object being set

168:   Output Parameter:
169: . pcomm - The ParallelComm for the DMMoab

171:   Level: beginner

173: @*/
174: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
175: {
177:   *pcomm = ((DM_Moab *)(dm)->data)->pcomm;
178:   return 0;
179: }

181: #endif /* MOAB_HAVE_MPI */

183: /*@C
184:   DMMoabSetInterface - Set the MOAB instance used with this DMMoab

186:   Collective

188:   Input Parameters:
189: + dm      - The DMMoab object being set
190: - mbiface - The MOAB instance being set on this DMMoab

192:   Level: beginner

194: @*/
195: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
196: {
197:   DM_Moab *dmmoab = (DM_Moab *)(dm)->data;

201: #ifdef MOAB_HAVE_MPI
202:   dmmoab->pcomm = NULL;
203: #endif
204:   dmmoab->mbiface          = mbiface;
205:   dmmoab->icreatedinstance = PETSC_FALSE;
206:   return 0;
207: }

209: /*@C
210:   DMMoabGetInterface - Get the MOAB instance used with this DMMoab

212:   Collective

214:   Input Parameter:
215: . dm      - The DMMoab object being set

217:   Output Parameter:
218: . mbiface - The MOAB instance set on this DMMoab

220:   Level: beginner

222: @*/
223: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
224: {
225:   static PetscBool cite = PETSC_FALSE;

228:   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, "
229:                                    "K. and Stimpson, C. and Ernst, C.},\n  year = {2004},  note = {Report}\n}\n",
230:                                    &cite));
231:   *mbiface = ((DM_Moab *)dm->data)->mbiface;
232:   return 0;
233: }

235: /*@C
236:   DMMoabSetLocalVertices - Set the entities having DOFs on this DMMoab

238:   Collective

240:   Input Parameters:
241: + dm    - The DMMoab object being set
242: - range - The entities treated by this DMMoab

244:   Level: beginner

246: @*/
247: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
248: {
249:   moab::Range tmpvtxs;
250:   DM_Moab    *dmmoab = (DM_Moab *)(dm)->data;

253:   dmmoab->vlocal->clear();
254:   dmmoab->vowned->clear();

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

258: #ifdef MOAB_HAVE_MPI
259:   moab::ErrorCode merr;
260:   /* filter based on parallel status */
261:   merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
262:   MBERRNM(merr);

264:   /* filter all the non-owned and shared entities out of the list */
265:   tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
266:   merr    = dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
267:   MBERRNM(merr);
268:   tmpvtxs         = moab::subtract(tmpvtxs, *dmmoab->vghost);
269:   *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
270: #else
271:   *dmmoab->vowned = *dmmoab->vlocal;
272: #endif

274:   /* compute and cache the sizes of local and ghosted entities */
275:   dmmoab->nloc   = dmmoab->vowned->size();
276:   dmmoab->nghost = dmmoab->vghost->size();
277: #ifdef MOAB_HAVE_MPI
278:   MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
279: #else
280:   dmmoab->n       = dmmoab->nloc;
281: #endif
282:   return 0;
283: }

285: /*@C
286:   DMMoabGetAllVertices - Get the entities having DOFs on this DMMoab

288:   Collective

290:   Input Parameter:
291: . dm    - The DMMoab object being set

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

296:   Level: beginner

298: @*/
299: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
300: {
302:   if (local) *local = *((DM_Moab *)dm->data)->vlocal;
303:   return 0;
304: }

306: /*@C
307:   DMMoabGetLocalVertices - Get the entities having DOFs on this DMMoab

309:   Collective

311:   Input Parameter:
312: . dm    - The DMMoab object being set

314:   Output Parameters:
315: + owned - The owned vertex entities in this DMMoab
316: - ghost - The ghosted entities (non-owned) stored locally in this partition

318:   Level: beginner

320: @*/
321: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
322: {
324:   if (owned) *owned = ((DM_Moab *)dm->data)->vowned;
325:   if (ghost) *ghost = ((DM_Moab *)dm->data)->vghost;
326:   return 0;
327: }

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

332:   Collective

334:   Input Parameter:
335: . dm    - The DMMoab object being set

337:   Output Parameter:
338: . range - The entities owned locally

340:   Level: beginner

342: @*/
343: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
344: {
346:   if (range) *range = ((DM_Moab *)dm->data)->elocal;
347:   return 0;
348: }

350: /*@C
351:   DMMoabSetLocalElements - Set the entities having DOFs on this DMMoab

353:   Collective

355:   Input Parameters:
356: + dm    - The DMMoab object being set
357: - range - The entities treated by this DMMoab

359:   Level: beginner

361: @*/
362: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
363: {
364:   DM_Moab *dmmoab = (DM_Moab *)(dm)->data;

367:   dmmoab->elocal->clear();
368:   dmmoab->eghost->clear();
369:   dmmoab->elocal->insert(range->begin(), range->end());
370: #ifdef MOAB_HAVE_MPI
371:   moab::ErrorCode merr;
372:   merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
373:   MBERRNM(merr);
374:   *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
375: #endif
376:   dmmoab->neleloc   = dmmoab->elocal->size();
377:   dmmoab->neleghost = dmmoab->eghost->size();
378: #ifdef MOAB_HAVE_MPI
379:   MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
380:   PetscInfo(dm, "Created %" PetscInt_FMT " local and %" PetscInt_FMT " global elements.\n", dmmoab->neleloc, dmmoab->nele);
381: #else
382:   dmmoab->nele    = dmmoab->neleloc;
383: #endif
384:   return 0;
385: }

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

390:   Collective

392:   Input Parameters:
393: + dm      - The DMMoab object being set
394: - ltogtag - The MOAB tag used for local to global ids

396:   Level: beginner

398: @*/
399: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
400: {
402:   ((DM_Moab *)dm->data)->ltog_tag = ltogtag;
403:   return 0;
404: }

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

409:   Collective

411:   Input Parameter:
412: . dm      - The DMMoab object being set

414:   Output Parameter:
415: . ltogtag - The MOAB tag used for local to global ids

417:   Level: beginner

419: @*/
420: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
421: {
423:   *ltog_tag = ((DM_Moab *)dm->data)->ltog_tag;
424:   return 0;
425: }

427: /*@C
428:   DMMoabSetBlockSize - Set the block size used with this DMMoab

430:   Collective

432:   Input Parameters:
433: + dm - The DMMoab object being set
434: - bs - The block size used with this DMMoab

436:   Level: beginner

438: @*/
439: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
440: {
442:   ((DM_Moab *)dm->data)->bs = bs;
443:   return 0;
444: }

446: /*@C
447:   DMMoabGetBlockSize - Get the block size used with this DMMoab

449:   Collective

451:   Input Parameter:
452: . dm - The DMMoab object being set

454:   Output Parameter:
455: . bs - The block size used with this DMMoab

457:   Level: beginner

459: @*/
460: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
461: {
463:   *bs = ((DM_Moab *)dm->data)->bs;
464:   return 0;
465: }

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

470:   Collective on dm

472:   Input Parameter:
473: . dm - The DMMoab object being set

475:   Output Parameters:
476: + neg - The number of global elements in the DMMoab instance
477: - nvg - The number of global vertices in the DMMoab instance

479:   Level: beginner

481: @*/
482: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
483: {
485:   if (neg) *neg = ((DM_Moab *)dm->data)->nele;
486:   if (nvg) *nvg = ((DM_Moab *)dm->data)->n;
487:   return 0;
488: }

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

493:   Collective on dm

495:   Input Parameter:
496: . dm - The DMMoab object being set

498:   Output Parameters:
499: + nel - The number of owned elements in this processor
500: . neg - The number of ghosted elements in this processor
501: . nvl - The number of owned vertices in this processor
502: - nvg - The number of ghosted vertices in this processor

504:   Level: beginner

506: @*/
507: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
508: {
510:   if (nel) *nel = ((DM_Moab *)dm->data)->neleloc;
511:   if (neg) *neg = ((DM_Moab *)dm->data)->neleghost;
512:   if (nvl) *nvl = ((DM_Moab *)dm->data)->nloc;
513:   if (nvg) *nvg = ((DM_Moab *)dm->data)->nghost;
514:   return 0;
515: }

517: /*@C
518:   DMMoabGetOffset - Get the local offset for the global vector

520:   Collective

522:   Input Parameter:
523: . dm - The DMMoab object being set

525:   Output Parameter:
526: . offset - The local offset for the global vector

528:   Level: beginner

530: @*/
531: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
532: {
534:   *offset = ((DM_Moab *)dm->data)->vstart;
535:   return 0;
536: }

538: /*@C
539:   DMMoabGetDimension - Get the dimension of the DM Mesh

541:   Collective

543:   Input Parameter:
544: . dm - The DMMoab object

546:   Output Parameter:
547: . dim - The dimension of DM

549:   Level: beginner

551: @*/
552: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
553: {
555:   *dim = ((DM_Moab *)dm->data)->dim;
556:   return 0;
557: }

559: /*@C
560:   DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
561:   generated through uniform refinement.

563:   Collective on dm

565:   Input Parameter:
566: . dm - The DMMoab object being set

568:   Output Parameter:
569: . nvg - The current mesh hierarchy level

571:   Level: beginner

573: @*/
574: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
575: {
577:   if (nlevel) *nlevel = ((DM_Moab *)dm->data)->hlevel;
578:   return 0;
579: }

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

584:   Collective

586:   Input Parameters:
587: + dm - The DMMoab object
588: - ehandle - The element entity handle

590:   Output Parameter:
591: . mat - The material ID for the current entity

593:   Level: beginner

595: @*/
596: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
597: {
598:   DM_Moab *dmmoab;

601:   if (*mat) {
602:     dmmoab = (DM_Moab *)(dm)->data;
603:     *mat   = dmmoab->materials[dmmoab->elocal->index(ehandle)];
604:   }
605:   return 0;
606: }

608: /*@C
609:   DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities

611:   Collective

613:   Input Parameters:
614: + dm - The DMMoab object
615: . nconn - Number of entities whose coordinates are needed
616: - conn - The vertex entity handles

618:   Output Parameter:
619: . vpos - The coordinates of the requested vertex entities

621:   Level: beginner

623: .seealso: `DMMoabGetVertexConnectivity()`
624: @*/
625: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
626: {
627:   DM_Moab        *dmmoab;
628:   moab::ErrorCode merr;

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

635:   /* Get connectivity information in MOAB canonical ordering */
636:   if (dmmoab->hlevel) {
637:     merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle *>(conn), nconn, dmmoab->hlevel, vpos);
638:     MBERRNM(merr);
639:   } else {
640:     merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);
641:     MBERRNM(merr);
642:   }
643:   return 0;
644: }

646: /*@C
647:   DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity

649:   Collective

651:   Input Parameters:
652: + dm - The DMMoab object
653: - vhandle - Vertex entity handle

655:   Output Parameters:
656: + nconn - Number of entities whose coordinates are needed
657: - conn - The vertex entity handles

659:   Level: beginner

661: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabRestoreVertexConnectivity()`
662: @*/
663: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt *nconn, moab::EntityHandle **conn)
664: {
665:   DM_Moab                        *dmmoab;
666:   std::vector<moab::EntityHandle> adj_entities, connect;
667:   moab::ErrorCode                 merr;

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

673:   /* Get connectivity information in MOAB canonical ordering */
674:   merr = dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION);
675:   MBERRNM(merr);
676:   merr = dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect);
677:   MBERRNM(merr);

679:   if (conn) {
680:     PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn);
681:     PetscArraycpy(*conn, &connect[0], connect.size());
682:   }
683:   if (nconn) *nconn = connect.size();
684:   return 0;
685: }

687: /*@C
688:   DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity

690:   Collective

692:   Input Parameters:
693: + dm - The DMMoab object
694: . vhandle - Vertex entity handle
695: . nconn - Number of entities whose coordinates are needed
696: - conn - The vertex entity handles

698:   Level: beginner

700: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`
701: @*/
702: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, moab::EntityHandle **conn)
703: {

707:   if (conn) PetscFree(*conn);
708:   if (nconn) *nconn = 0;
709:   return 0;
710: }

712: /*@C
713:   DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity

715:   Collective

717:   Input Parameters:
718: + dm - The DMMoab object
719: - ehandle - Vertex entity handle

721:   Output Parameters:
722: + nconn - Number of entities whose coordinates are needed
723: - conn - The vertex entity handles

725:   Level: beginner

727: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`, `DMMoabRestoreVertexConnectivity()`
728: @*/
729: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, const moab::EntityHandle **conn)
730: {
731:   DM_Moab                        *dmmoab;
732:   const moab::EntityHandle       *connect;
733:   std::vector<moab::EntityHandle> vconn;
734:   moab::ErrorCode                 merr;
735:   PetscInt                        nnodes;

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

741:   /* Get connectivity information in MOAB canonical ordering */
742:   merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes);
743:   MBERRNM(merr);
744:   if (conn) *conn = connect;
745:   if (nconn) *nconn = nnodes;
746:   return 0;
747: }

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

752:   Collective

754:   Input Parameters:
755: + dm - The DMMoab object
756: - ent - Entity handle

758:   Output Parameter:
759: . ent_on_boundary - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise

761:   Level: beginner

763: .seealso: `DMMoabCheckBoundaryVertices()`
764: @*/
765: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool *ent_on_boundary)
766: {
767:   moab::EntityType etype;
768:   DM_Moab         *dmmoab;
769:   PetscInt         edim;

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

775:   /* get the entity type and handle accordingly */
776:   etype = dmmoab->mbiface->type_from_handle(ent);

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

782:   *ent_on_boundary = PETSC_FALSE;
783:   if (etype == moab::MBVERTEX && edim == 0) {
784:     *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
785:   } else {
786:     if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
787:       if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
788:     } else { /* next check the lower-dimensional faces */
789:       if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
790:     }
791:   }
792:   return 0;
793: }

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

798:   Input Parameters:
799: + dm - The DMMoab object
800: . nconn - Number of handles
801: - cnt - Array of entity handles

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

806:   Level: beginner

808: .seealso: `DMMoabIsEntityOnBoundary()`
809: @*/
810: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool *isbdvtx)
811: {
812:   DM_Moab *dmmoab;
813:   PetscInt i;

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

820:   for (i = 0; i < nconn; ++i) isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
821:   return 0;
822: }

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

827:   Input Parameter:
828: . dm - The DMMoab object

830:   Output Parameters:
831: + bdvtx - Boundary vertices
832: . bdelems - Boundary elements
833: - bdfaces - Boundary faces

835:   Level: beginner

837: .seealso: `DMMoabCheckBoundaryVertices()`, `DMMoabIsEntityOnBoundary()`
838: @*/
839: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range **bdelems, const moab::Range **bdfaces)
840: {
841:   DM_Moab *dmmoab;

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

846:   if (bdvtx) *bdvtx = dmmoab->bndyvtx;
847:   if (bdfaces) *bdfaces = dmmoab->bndyfaces;
848:   if (bdelems) *bdfaces = dmmoab->bndyelems;
849:   return 0;
850: }

852: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
853: {
854:   PetscInt        i;
855:   moab::ErrorCode merr;
856:   DM_Moab        *dmmoab = (DM_Moab *)dm->data;


860:   dmmoab->refct--;
861:   if (!dmmoab->refct) {
862:     delete dmmoab->vlocal;
863:     delete dmmoab->vowned;
864:     delete dmmoab->vghost;
865:     delete dmmoab->elocal;
866:     delete dmmoab->eghost;
867:     delete dmmoab->bndyvtx;
868:     delete dmmoab->bndyfaces;
869:     delete dmmoab->bndyelems;

871:     PetscFree(dmmoab->gsindices);
872:     PetscFree2(dmmoab->gidmap, dmmoab->lidmap);
873:     PetscFree(dmmoab->dfill);
874:     PetscFree(dmmoab->ofill);
875:     PetscFree(dmmoab->materials);
876:     if (dmmoab->fieldNames) {
877:       for (i = 0; i < dmmoab->numFields; i++) PetscFree(dmmoab->fieldNames[i]);
878:       PetscFree(dmmoab->fieldNames);
879:     }

881:     if (dmmoab->nhlevels) {
882:       PetscFree(dmmoab->hsets);
883:       dmmoab->nhlevels = 0;
884:       if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
885:       dmmoab->hierarchy = NULL;
886:     }

888:     if (dmmoab->icreatedinstance) {
889:       delete dmmoab->pcomm;
890:       merr = dmmoab->mbiface->delete_mesh();
891:       MBERRNM(merr);
892:       delete dmmoab->mbiface;
893:     }
894:     dmmoab->mbiface = NULL;
895: #ifdef MOAB_HAVE_MPI
896:     dmmoab->pcomm = NULL;
897: #endif
898:     VecScatterDestroy(&dmmoab->ltog_sendrecv);
899:     ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map);
900:     PetscFree(dm->data);
901:   }
902:   return 0;
903: }

905: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(DM dm, PetscOptionItems *PetscOptionsObject)
906: {
907:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

909:   PetscOptionsHeadBegin(PetscOptionsObject, "DMMoab Options");
910:   PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL, 0);
911:   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);
912:   /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
913:   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);
914:   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);
915:   PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum *)&dmmoab->read_mode, NULL);
916:   PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum *)&dmmoab->write_mode, NULL);
917:   PetscOptionsHeadEnd();
918:   return 0;
919: }

921: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
922: {
923:   moab::ErrorCode       merr;
924:   Vec                   local, global;
925:   IS                    from, to;
926:   moab::Range::iterator iter;
927:   PetscInt              i, j, f, bs, vent, totsize, *lgmap;
928:   DM_Moab              *dmmoab = (DM_Moab *)dm->data;
929:   moab::Range           adjs;

932:   /* Get the local and shared vertices and cache it */
934: #ifdef MOAB_HAVE_MPI
936: #endif

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

944: #ifdef MOAB_HAVE_MPI
945:     /* filter based on parallel status */
946:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
947:     MBERRNM(merr);

949:     /* filter all the non-owned and shared entities out of the list */
950:     // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
951:     adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
952:     merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
953:     MBERRNM(merr);
954:     adjs            = moab::subtract(adjs, *dmmoab->vghost);
955:     *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
956: #else
957:     *dmmoab->vowned = *dmmoab->vlocal;
958: #endif

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

964: #ifdef MOAB_HAVE_MPI
965:     MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
966:     PetscInfo(NULL, "Filset ID: %lu, Vertices: local - %zu, owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->fileset, dmmoab->vlocal->size(), dmmoab->nloc, dmmoab->nghost);
967: #else
968:     dmmoab->n       = dmmoab->nloc;
969: #endif
970:   }

972:   {
973:     /* get the information about the local elements in the mesh */
974:     dmmoab->eghost->clear();

976:     /* first decipher the leading dimension */
977:     for (i = 3; i > 0; i--) {
978:       dmmoab->elocal->clear();
979:       merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false);
980:       MBERRNM(merr);

982:       /* store the current mesh dimension */
983:       if (dmmoab->elocal->size()) {
984:         dmmoab->dim = i;
985:         break;
986:       }
987:     }

989:     DMSetDimension(dm, dmmoab->dim);

991: #ifdef MOAB_HAVE_MPI
992:     /* filter the ghosted and owned element list */
993:     *dmmoab->eghost = *dmmoab->elocal;
994:     merr            = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
995:     MBERRNM(merr);
996:     *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
997: #endif

999:     dmmoab->neleloc   = dmmoab->elocal->size();
1000:     dmmoab->neleghost = dmmoab->eghost->size();

1002: #ifdef MOAB_HAVE_MPI
1003:     MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
1004:     PetscInfo(NULL, "%d-dim elements: owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost);
1005: #else
1006:     dmmoab->nele    = dmmoab->neleloc;
1007: #endif
1008:   }

1010:   bs = dmmoab->bs;
1011:   if (!dmmoab->ltog_tag) {
1012:     /* Get the global ID tag. The global ID tag is applied to each
1013:        vertex. It acts as an global identifier which MOAB uses to
1014:        assemble the individual pieces of the mesh */
1015:     merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
1016:     MBERRNM(merr);
1017:   }

1019:   totsize = dmmoab->vlocal->size();
1021:   PetscCalloc1(totsize, &dmmoab->gsindices);
1022:   {
1023:     /* first get the local indices */
1024:     merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]);
1025:     MBERRNM(merr);
1026:     if (dmmoab->nghost) { /* next get the ghosted indices */
1027:       merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]);
1028:       MBERRNM(merr);
1029:     }

1031:     /* find out the local and global minima of GLOBAL_ID */
1032:     dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1033:     for (i = 0; i < totsize; ++i) {
1034:       if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1035:       if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1036:     }

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

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

1046:     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]);
1047:   }
1049:              dmmoab->numFields);

1051:   {
1052:     dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1053:     dmmoab->seqend   = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1054:     PetscInfo(NULL, "SEQUENCE: Local [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "]\n", dmmoab->seqstart, dmmoab->seqend);

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

1059:     i = j = 0;
1060:     /* set the owned vertex data first */
1061:     for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1062:       vent                 = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1063:       dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1064:       dmmoab->lidmap[vent] = i;
1065:       for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1066:     }
1067:     /* next arrange all the ghosted data information */
1068:     for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1069:       vent                 = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1070:       dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1071:       dmmoab->lidmap[vent] = i;
1072:       for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1073:     }

1075:     /* We need to create the Global to Local Vector Scatter Contexts
1076:        1) First create a local and global vector
1077:        2) Create a local and global IS
1078:        3) Create VecScatter and LtoGMapping objects
1079:        4) Cleanup the IS and Vec objects
1080:     */
1081:     DMCreateGlobalVector(dm, &global);
1082:     DMCreateLocalVector(dm, &local);

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

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

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

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

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

1101:     /* clean up IS, Vec */
1102:     PetscFree(lgmap);
1103:     ISDestroy(&from);
1104:     ISDestroy(&to);
1105:     VecDestroy(&local);
1106:     VecDestroy(&global);
1107:   }

1109:   dmmoab->bndyvtx   = new moab::Range();
1110:   dmmoab->bndyfaces = new moab::Range();
1111:   dmmoab->bndyelems = new moab::Range();
1112:   /* skin the boundary and store nodes */
1113:   if (!dmmoab->hlevel) {
1114:     /* get the skin vertices of boundary faces for the current partition and then filter
1115:        the local, boundary faces, vertices and elements alone via PSTATUS flags;
1116:        this should not give us any ghosted boundary, but if user needs such a functionality
1117:        it would be easy to add it based on the find_skin query below */
1118:     moab::Skinner skinner(dmmoab->mbiface);

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

1124: #ifdef MOAB_HAVE_MPI
1125:     /* filter all the non-owned and shared entities out of the list */
1126:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1127:     MBERRNM(merr);
1128:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT);
1129:     MBERRNM(merr);
1130: #endif

1132:     /* get all the nodes via connectivity and the parent elements via adjacency information */
1133:     merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false);
1134:     MBERRNM(merr);
1135:     merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION);
1136:     MBERRNM(merr);
1137:   } else {
1138:     /* Let us query the hierarchy manager and get the results directly for this level */
1139:     for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1140:       moab::EntityHandle elemHandle = *iter;
1141:       if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1142:         dmmoab->bndyelems->insert(elemHandle);
1143:         /* For this boundary element, query the vertices and add them to the list */
1144:         std::vector<moab::EntityHandle> connect;
1145:         merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect);
1146:         MBERRNM(merr);
1147:         for (unsigned iv = 0; iv < connect.size(); ++iv)
1148:           if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv])) dmmoab->bndyvtx->insert(connect[iv]);
1149:         /* Next, let us query the boundary faces and add them also to the list */
1150:         std::vector<moab::EntityHandle> faces;
1151:         merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim - 1, faces);
1152:         MBERRNM(merr);
1153:         for (unsigned ifa = 0; ifa < faces.size(); ++ifa)
1154:           if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa])) dmmoab->bndyfaces->insert(faces[ifa]);
1155:       }
1156:     }
1157: #ifdef MOAB_HAVE_MPI
1158:     /* filter all the non-owned and shared entities out of the list */
1159:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1160:     MBERRNM(merr);
1161:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1162:     MBERRNM(merr);
1163:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1164:     MBERRNM(merr);
1165: #endif
1166:   }
1167:   PetscInfo(NULL, "Found %zu boundary vertices, %zu boundary faces and %zu boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size());

1169:   /* Get the material sets and populate the data for all locally owned elements */
1170:   {
1171:     PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials);
1172:     /* Get the count of entities of particular type from dmmoab->elocal
1173:        -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1174:     moab::Range msets;
1175:     merr = dmmoab->mbiface->get_entities_by_type_and_tag(dmmoab->fileset, moab::MBENTITYSET, &dmmoab->material_tag, NULL, 1, msets, moab::Interface::UNION);
1176:     MB_CHK_ERR(merr);
1177:     if (msets.size() == 0) PetscInfo(NULL, "No material sets found in the fileset.");

1179:     for (unsigned i = 0; i < msets.size(); ++i) {
1180:       moab::Range msetelems;
1181:       merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);
1182:       MB_CHK_ERR(merr);
1183: #ifdef MOAB_HAVE_MPI
1184:       /* filter all the non-owned and shared entities out of the list */
1185:       merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1186:       MBERRNM(merr);
1187: #endif

1189:       int                partID;
1190:       moab::EntityHandle mset = msets[i];
1191:       merr                    = dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID);
1192:       MB_CHK_ERR(merr);

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

1198:   return 0;
1199: }

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

1204:   Collective

1206:   Input Parameters:
1207: + dm - The DM object
1208: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1209: . conn - The connectivity of the element
1210: - nverts - The number of vertices that form the element

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

1215:   Level: beginner

1217: .seealso: `DMMoabCreateSubmesh()`, `DMMoabCreateElement()`
1218: @*/
1219: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal *coords, PetscInt nverts, moab::Range *overts)
1220: {
1221:   moab::ErrorCode merr;
1222:   DM_Moab        *dmmoab;
1223:   moab::Range     verts;


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

1230:   /* Insert new points */
1231:   merr = dmmoab->mbiface->create_vertices(&coords[0], nverts, verts);
1232:   MBERRNM(merr);
1233:   merr = dmmoab->mbiface->add_entities(dmmoab->fileset, verts);
1234:   MBERRNM(merr);

1236:   if (overts) *overts = verts;
1237:   return 0;
1238: }

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

1243:   Collective

1245:   Input Parameters:
1246: + dm - The DM object
1247: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1248: . conn - The connectivity of the element
1249: - nverts - The number of vertices that form the element

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

1254:   Level: beginner

1256: .seealso: `DMMoabCreateSubmesh()`, `DMMoabCreateVertices()`
1257: @*/
1258: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle *conn, PetscInt nverts, moab::EntityHandle *oelem)
1259: {
1260:   moab::ErrorCode    merr;
1261:   DM_Moab           *dmmoab;
1262:   moab::EntityHandle elem;


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

1269:   /* Insert new element */
1270:   merr = dmmoab->mbiface->create_element(type, conn, nverts, elem);
1271:   MBERRNM(merr);
1272:   merr = dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1);
1273:   MBERRNM(merr);

1275:   if (oelem) *oelem = elem;
1276:   return 0;
1277: }

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

1284:   Collective

1286:   Input Parameters:
1287: . dm - The DM object

1289:   Output Parameter:
1290: . newdm  - The sub DM object with updated set information

1292:   Level: advanced

1294: .seealso: `DMCreate()`, `DMMoabCreateVertices()`, `DMMoabCreateElement()`
1295: @*/
1296: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1297: {
1298:   DM_Moab        *dmmoab;
1299:   DM_Moab        *ndmmoab;
1300:   moab::ErrorCode merr;


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

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

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

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

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

1319:   /* create a meshset and then add old fileset as child */
1320:   merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal);
1321:   MBERR("Adding child vertices to parent failed", merr);
1322:   merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal);
1323:   MBERR("Adding child elements to parent failed", merr);

1325:   /* preserve the field association between the parent and sub-mesh objects */
1326:   DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames);
1327:   return 0;
1328: }

1330: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1331: {
1332:   DM_Moab    *dmmoab = (DM_Moab *)(dm)->data;
1333:   const char *name;
1334:   MPI_Comm    comm;
1335:   PetscMPIInt size;

1337:   PetscObjectGetComm((PetscObject)dm, &comm);
1338:   MPI_Comm_size(comm, &size);
1339:   PetscObjectGetName((PetscObject)dm, &name);
1340:   PetscViewerASCIIPushTab(viewer);
1341:   if (name) PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimensions:\n", name, dmmoab->dim);
1342:   else PetscViewerASCIIPrintf(viewer, "Mesh in %" PetscInt_FMT " dimensions:\n", dmmoab->dim);
1343:   /* print details about the global mesh */
1344:   {
1345:     PetscViewerASCIIPushTab(viewer);
1346:     PetscViewerASCIIPrintf(viewer, "Sizes: cells=%" PetscInt_FMT ", vertices=%" PetscInt_FMT ", blocks=%" PetscInt_FMT "\n", dmmoab->nele, dmmoab->n, dmmoab->bs);
1347:     /* print boundary data */
1348:     PetscViewerASCIIPrintf(viewer, "Boundary trace:\n");
1349:     {
1350:       PetscViewerASCIIPushTab(viewer);
1351:       PetscViewerASCIIPrintf(viewer, "cells=%zu, faces=%zu, vertices=%zu\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size());
1352:       PetscViewerASCIIPopTab(viewer);
1353:     }
1354:     /* print field data */
1355:     PetscViewerASCIIPrintf(viewer, "Fields: %" PetscInt_FMT " components\n", dmmoab->numFields);
1356:     {
1357:       PetscViewerASCIIPushTab(viewer);
1358:       for (int i = 0; i < dmmoab->numFields; ++i) PetscViewerASCIIPrintf(viewer, "[%" PetscInt_FMT "] - %s\n", i, dmmoab->fieldNames[i]);
1359:       PetscViewerASCIIPopTab(viewer);
1360:     }
1361:     PetscViewerASCIIPopTab(viewer);
1362:   }
1363:   PetscViewerASCIIPopTab(viewer);
1364:   PetscViewerFlush(viewer);
1365:   return 0;
1366: }

1368: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1369: {
1370:   return 0;
1371: }

1373: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1374: {
1375:   return 0;
1376: }

1378: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1379: {
1380:   PetscBool iascii, ishdf5, isvtk;

1384:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1385:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);
1386:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
1387:   if (iascii) {
1388:     DMMoabView_Ascii(dm, viewer);
1389:   } else if (ishdf5) {
1390: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1391:     PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
1392:     DMMoabView_HDF5(dm, viewer);
1393:     PetscViewerPopFormat(viewer);
1394: #else
1395:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1396: #endif
1397:   } else if (isvtk) {
1398:     DMMoabView_VTK(dm, viewer);
1399:   }
1400:   return 0;
1401: }

1403: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
1404: {
1405:   dm->ops->view                     = DMView_Moab;
1406:   dm->ops->load                     = NULL /* DMLoad_Moab */;
1407:   dm->ops->setfromoptions           = DMSetFromOptions_Moab;
1408:   dm->ops->clone                    = DMClone_Moab;
1409:   dm->ops->setup                    = DMSetUp_Moab;
1410:   dm->ops->createlocalsection       = NULL;
1411:   dm->ops->createdefaultconstraints = NULL;
1412:   dm->ops->createglobalvector       = DMCreateGlobalVector_Moab;
1413:   dm->ops->createlocalvector        = DMCreateLocalVector_Moab;
1414:   dm->ops->getlocaltoglobalmapping  = NULL;
1415:   dm->ops->createfieldis            = NULL;
1416:   dm->ops->createcoordinatedm       = NULL /* DMCreateCoordinateDM_Moab */;
1417:   dm->ops->getcoloring              = NULL;
1418:   dm->ops->creatematrix             = DMCreateMatrix_Moab;
1419:   dm->ops->createinterpolation      = DMCreateInterpolation_Moab;
1420:   dm->ops->createinjection          = NULL /* DMCreateInjection_Moab */;
1421:   dm->ops->refine                   = DMRefine_Moab;
1422:   dm->ops->coarsen                  = DMCoarsen_Moab;
1423:   dm->ops->refinehierarchy          = DMRefineHierarchy_Moab;
1424:   dm->ops->coarsenhierarchy         = DMCoarsenHierarchy_Moab;
1425:   dm->ops->globaltolocalbegin       = DMGlobalToLocalBegin_Moab;
1426:   dm->ops->globaltolocalend         = DMGlobalToLocalEnd_Moab;
1427:   dm->ops->localtoglobalbegin       = DMLocalToGlobalBegin_Moab;
1428:   dm->ops->localtoglobalend         = DMLocalToGlobalEnd_Moab;
1429:   dm->ops->destroy                  = DMDestroy_Moab;
1430:   dm->ops->createsubdm              = NULL /* DMCreateSubDM_Moab */;
1431:   dm->ops->getdimpoints             = NULL /* DMGetDimPoints_Moab */;
1432:   dm->ops->locatepoints             = NULL /* DMLocatePoints_Moab */;
1433:   return 0;
1434: }

1436: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1437: {
1438:   /* get all the necessary handles from the private DM object */
1439:   (*newdm)->data = (DM_Moab *)dm->data;
1440:   ((DM_Moab *)dm->data)->refct++;

1442:   PetscObjectChangeTypeName((PetscObject)*newdm, DMMOAB);
1443:   DMInitialize_Moab(*newdm);
1444:   return 0;
1445: }

1447: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1448: {
1450:   PetscNew((DM_Moab **)&dm->data);

1452:   ((DM_Moab *)dm->data)->bs            = 1;
1453:   ((DM_Moab *)dm->data)->numFields     = 1;
1454:   ((DM_Moab *)dm->data)->n             = 0;
1455:   ((DM_Moab *)dm->data)->nloc          = 0;
1456:   ((DM_Moab *)dm->data)->nghost        = 0;
1457:   ((DM_Moab *)dm->data)->nele          = 0;
1458:   ((DM_Moab *)dm->data)->neleloc       = 0;
1459:   ((DM_Moab *)dm->data)->neleghost     = 0;
1460:   ((DM_Moab *)dm->data)->ltog_map      = NULL;
1461:   ((DM_Moab *)dm->data)->ltog_sendrecv = NULL;

1463:   ((DM_Moab *)dm->data)->refct  = 1;
1464:   ((DM_Moab *)dm->data)->parent = NULL;
1465:   ((DM_Moab *)dm->data)->vlocal = new moab::Range();
1466:   ((DM_Moab *)dm->data)->vowned = new moab::Range();
1467:   ((DM_Moab *)dm->data)->vghost = new moab::Range();
1468:   ((DM_Moab *)dm->data)->elocal = new moab::Range();
1469:   ((DM_Moab *)dm->data)->eghost = new moab::Range();

1471:   DMInitialize_Moab(dm);
1472:   return 0;
1473: }