Actual source code: dmmoab.cxx

petsc-master 2019-08-19
Report Typos and Errors
  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);


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

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

 56:   Collective

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

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

 64:   Level: beginner

 66: @*/
 67: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
 68: {

 73:   DMCreate(comm, dmb);
 74:   DMSetType(*dmb, DMMOAB);
 75:   return(0);
 76: }

 78: /*@
 79:   DMMoabCreateMoab - Creates a DMMoab object, optionally from an instance and other data

 81:   Collective

 83:   Input Parameter:
 84: + comm - The communicator for the DMMoab object
 85: . mbiface - (ptr to) the MOAB Instance; if passed in NULL, MOAB instance is created inside PETSc, and destroyed
 86:          along with the DMMoab
 87: . pcomm - (ptr to) a ParallelComm; if NULL, creates one internally for the whole communicator
 88: . ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
 89: - range - If non-NULL, contains range of entities to which DOFs will be assigned

 91:   Output Parameter:
 92: . dmb  - The DMMoab object

 94:   Level: intermediate

 96: @*/
 97: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
 98: {
100:   moab::ErrorCode merr;
101:   DM             dmmb;
102:   DM_Moab        *dmmoab;


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

110:   if (!mbiface) {
111:     dmmoab->mbiface = new moab::Core();
112:     dmmoab->icreatedinstance = PETSC_TRUE;
113:   }
114:   else {
115:     dmmoab->mbiface = mbiface;
116:     dmmoab->icreatedinstance = PETSC_FALSE;
117:   }

119:   /* by default the fileset = root set. This set stores the hierarchy of entities belonging to current DM */
120:   dmmoab->fileset = 0;
121:   dmmoab->hlevel = 0;
122:   dmmoab->nghostrings = 0;

124: #ifdef MOAB_HAVE_MPI
125:   moab::EntityHandle partnset;

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

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

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

147:   /* set global ID tag handle */
148:   if (ltog_tag && *ltog_tag) {
149:     DMMoabSetLocalToGlobalTag(dmmb, *ltog_tag);
150:   }
151:   else {
152:     merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag); MBERRNM(merr);
153:     if (ltog_tag) *ltog_tag = dmmoab->ltog_tag;
154:   }

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

158:   /* set the local range of entities (vertices) of interest */
159:   if (range) {
160:     DMMoabSetLocalVertices(dmmb, range);
161:   }
162:   *dmb = dmmb;
163:   return(0);
164: }


167: #ifdef MOAB_HAVE_MPI

169: /*@
170:   DMMoabGetParallelComm - Get the ParallelComm used with this DMMoab

172:   Collective

174:   Input Parameter:
175: . dm    - The DMMoab object being set

177:   Output Parameter:
178: . pcomm - The ParallelComm for the DMMoab

180:   Level: beginner

182: @*/
183: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
184: {
187:   *pcomm = ((DM_Moab*)(dm)->data)->pcomm;
188:   return(0);
189: }

191: #endif /* MOAB_HAVE_MPI */


194: /*@
195:   DMMoabSetInterface - Set the MOAB instance used with this DMMoab

197:   Collective

199:   Input Parameter:
200: + dm      - The DMMoab object being set
201: - mbiface - The MOAB instance being set on this DMMoab

203:   Level: beginner

205: @*/
206: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
207: {
208:   DM_Moab        *dmmoab = (DM_Moab*)(dm)->data;

213: #ifdef MOAB_HAVE_MPI
214:   dmmoab->pcomm = NULL;
215: #endif
216:   dmmoab->mbiface = mbiface;
217:   dmmoab->icreatedinstance = PETSC_FALSE;
218:   return(0);
219: }


222: /*@
223:   DMMoabGetInterface - Get the MOAB instance used with this DMMoab

225:   Collective

227:   Input Parameter:
228: . dm      - The DMMoab object being set

230:   Output Parameter:
231: . mbiface - The MOAB instance set on this DMMoab

233:   Level: beginner

235: @*/
236: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
237: {
238:   PetscErrorCode   ierr;
239:   static PetscBool cite = PETSC_FALSE;

243:   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, K. and Stimpson, C. and Ernst, C.},\n  year = {2004},  note = {Report}\n}\n", &cite);
244:   *mbiface = ((DM_Moab*)dm->data)->mbiface;
245:   return(0);
246: }


249: /*@
250:   DMMoabSetLocalVertices - Set the entities having DOFs on this DMMoab

252:   Collective

254:   Input Parameter:
255: + dm    - The DMMoab object being set
256: - range - The entities treated by this DMMoab

258:   Level: beginner

260: @*/
261: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
262: {
263:   moab::Range     tmpvtxs;
264:   DM_Moab        *dmmoab = (DM_Moab*)(dm)->data;

268:   dmmoab->vlocal->clear();
269:   dmmoab->vowned->clear();

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

273: #ifdef MOAB_HAVE_MPI
274:   moab::ErrorCode merr;
275:   /* filter based on parallel status */
276:   merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned); MBERRNM(merr);

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

287:   /* compute and cache the sizes of local and ghosted entities */
288:   dmmoab->nloc = dmmoab->vowned->size();
289:   dmmoab->nghost = dmmoab->vghost->size();
290: #ifdef MOAB_HAVE_MPI
291:   PetscErrorCode  ierr;
292:   MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
293: #else
294:   dmmoab->n = dmmoab->nloc;
295: #endif
296:   return(0);
297: }


300: /*@
301:   DMMoabGetAllVertices - Get the entities having DOFs on this DMMoab

303:   Collective

305:   Input Parameter:
306: . dm    - The DMMoab object being set

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

311:   Level: beginner

313: @*/
314: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
315: {
318:   if (local) *local = *((DM_Moab*)dm->data)->vlocal;
319:   return(0);
320: }



324: /*@
325:   DMMoabGetLocalVertices - Get the entities having DOFs on this DMMoab

327:   Collective

329:   Input Parameter:
330: . dm    - The DMMoab object being set

332:   Output Parameters:
333: + owned - The owned vertex entities in this DMMoab
334: - ghost - The ghosted entities (non-owned) stored locally in this partition

336:   Level: beginner

338: @*/
339: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
340: {
343:   if (owned) *owned = ((DM_Moab*)dm->data)->vowned;
344:   if (ghost) *ghost = ((DM_Moab*)dm->data)->vghost;
345:   return(0);
346: }

348: /*@
349:   DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned

351:   Collective

353:   Input Parameter:
354: . dm    - The DMMoab object being set

356:   Output Parameter:
357: . range - The entities owned locally

359:   Level: beginner

361: @*/
362: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
363: {
366:   if (range) *range = ((DM_Moab*)dm->data)->elocal;
367:   return(0);
368: }


371: /*@
372:   DMMoabSetLocalElements - Set the entities having DOFs on this DMMoab

374:   Collective

376:   Input Parameters:
377: + dm    - The DMMoab object being set
378: - range - The entities treated by this DMMoab

380:   Level: beginner

382: @*/
383: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
384: {
385:   DM_Moab        *dmmoab = (DM_Moab*)(dm)->data;

389:   dmmoab->elocal->clear();
390:   dmmoab->eghost->clear();
391:   dmmoab->elocal->insert(range->begin(), range->end());
392: #ifdef MOAB_HAVE_MPI
393:   moab::ErrorCode merr;
394:   merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
395:   *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
396: #endif
397:   dmmoab->neleloc = dmmoab->elocal->size();
398:   dmmoab->neleghost = dmmoab->eghost->size();
399: #ifdef MOAB_HAVE_MPI
400:   PetscErrorCode  ierr;
401:   MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
402:   PetscInfo2(dm, "Created %D local and %D global elements.\n", dmmoab->neleloc, dmmoab->nele);
403: #else
404:   dmmoab->nele = dmmoab->neleloc;
405: #endif
406:   return(0);
407: }


410: /*@
411:   DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering

413:   Collective

415:   Input Parameters:
416: + dm      - The DMMoab object being set
417: - ltogtag - The MOAB tag used for local to global ids

419:   Level: beginner

421: @*/
422: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
423: {
426:   ((DM_Moab*)dm->data)->ltog_tag = ltogtag;
427:   return(0);
428: }


431: /*@
432:   DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering

434:   Collective

436:   Input Parameter:
437: . dm      - The DMMoab object being set

439:   Output Parameter:
440: . ltogtag - The MOAB tag used for local to global ids

442:   Level: beginner

444: @*/
445: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
446: {
449:   *ltog_tag = ((DM_Moab*)dm->data)->ltog_tag;
450:   return(0);
451: }


454: /*@
455:   DMMoabSetBlockSize - Set the block size used with this DMMoab

457:   Collective

459:   Input Parameter:
460: + dm - The DMMoab object being set
461: - bs - The block size used with this DMMoab

463:   Level: beginner

465: @*/
466: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
467: {
470:   ((DM_Moab*)dm->data)->bs = bs;
471:   return(0);
472: }


475: /*@
476:   DMMoabGetBlockSize - Get the block size used with this DMMoab

478:   Collective

480:   Input Parameter:
481: . dm - The DMMoab object being set

483:   Output Parameter:
484: . bs - The block size used with this DMMoab

486:   Level: beginner

488: @*/
489: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
490: {
493:   *bs = ((DM_Moab*)dm->data)->bs;
494:   return(0);
495: }


498: /*@
499:   DMMoabGetSize - Get the global vertex size used with this DMMoab

501:   Collective on dm

503:   Input Parameter:
504: . dm - The DMMoab object being set

506:   Output Parameter:
507: + neg - The number of global elements in the DMMoab instance
508: - nvg - The number of global vertices in the DMMoab instance

510:   Level: beginner

512: @*/
513: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
514: {
517:   if (neg) *neg = ((DM_Moab*)dm->data)->nele;
518:   if (nvg) *nvg = ((DM_Moab*)dm->data)->n;
519:   return(0);
520: }


523: /*@
524:   DMMoabGetLocalSize - Get the local and ghosted vertex size used with this DMMoab

526:   Collective on dm

528:   Input Parameter:
529: . dm - The DMMoab object being set

531:   Output Parameter:
532: + nel - The number of owned elements in this processor
533: . neg - The number of ghosted elements in this processor
534: . nvl - The number of owned vertices in this processor
535: - nvg - The number of ghosted vertices in this processor

537:   Level: beginner

539: @*/
540: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
541: {
544:   if (nel) *nel = ((DM_Moab*)dm->data)->neleloc;
545:   if (neg) *neg = ((DM_Moab*)dm->data)->neleghost;
546:   if (nvl) *nvl = ((DM_Moab*)dm->data)->nloc;
547:   if (nvg) *nvg = ((DM_Moab*)dm->data)->nghost;
548:   return(0);
549: }


552: /*@
553:   DMMoabGetOffset - Get the local offset for the global vector

555:   Collective

557:   Input Parameter:
558: . dm - The DMMoab object being set

560:   Output Parameter:
561: . offset - The local offset for the global vector

563:   Level: beginner

565: @*/
566: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
567: {
570:   *offset = ((DM_Moab*)dm->data)->vstart;
571:   return(0);
572: }


575: /*@
576:   DMMoabGetDimension - Get the dimension of the DM Mesh

578:   Collective

580:   Input Parameter:
581: . dm - The DMMoab object

583:   Output Parameter:
584: . dim - The dimension of DM

586:   Level: beginner

588: @*/
589: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
590: {
593:   *dim = ((DM_Moab*)dm->data)->dim;
594:   return(0);
595: }


598: /*@
599:   DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
600:   generated through uniform refinement.

602:   Collective on dm

604:   Input Parameter:
605: . dm - The DMMoab object being set

607:   Output Parameter:
608: . nvg - The current mesh hierarchy level

610:   Level: beginner

612: @*/
613: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
614: {
617:   if (nlevel) *nlevel = ((DM_Moab*)dm->data)->hlevel;
618:   return(0);
619: }


622: /*@
623:   DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh

625:   Collective

627:   Input Parameter:
628: + dm - The DMMoab object
629: - ehandle - The element entity handle

631:   Output Parameter:
632: . mat - The material ID for the current entity

634:   Level: beginner

636: @*/
637: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
638: {
639:   DM_Moab         *dmmoab;

643:   if (*mat) {
644:     dmmoab = (DM_Moab*)(dm)->data;
645:     *mat = dmmoab->materials[dmmoab->elocal->index(ehandle)];
646:   }
647:   return(0);
648: }


651: /*@
652:   DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities

654:   Collective

656:   Input Parameter:
657: + dm - The DMMoab object
658: . nconn - Number of entities whose coordinates are needed
659: - conn - The vertex entity handles

661:   Output Parameter:
662: . vpos - The coordinates of the requested vertex entities

664:   Level: beginner

666: .seealso: DMMoabGetVertexConnectivity()
667: @*/
668: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
669: {
670:   DM_Moab         *dmmoab;
671:   moab::ErrorCode merr;

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

679:   /* Get connectivity information in MOAB canonical ordering */
680:   if (dmmoab->hlevel) {
681:     merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle*>(conn), nconn, dmmoab->hlevel, vpos);MBERRNM(merr);
682:   }
683:   else {
684:     merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);MBERRNM(merr);
685:   }
686:   return(0);
687: }


690: /*@
691:   DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity

693:   Collective

695:   Input Parameter:
696: + dm - The DMMoab object
697: - vhandle - Vertex entity handle

699:   Output Parameter:
700: + nconn - Number of entities whose coordinates are needed
701: - conn - The vertex entity handles

703:   Level: beginner

705: .seealso: DMMoabGetVertexCoordinates(), DMMoabRestoreVertexConnectivity()
706: @*/
707: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt* nconn, moab::EntityHandle **conn)
708: {
709:   DM_Moab        *dmmoab;
710:   std::vector<moab::EntityHandle> adj_entities, connect;
711:   PetscErrorCode  ierr;
712:   moab::ErrorCode merr;

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

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

723:   if (conn) {
724:     PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn);
725:     PetscArraycpy(*conn, &connect[0], connect.size());
726:   }
727:   if (nconn) *nconn = connect.size();
728:   return(0);
729: }


732: /*@
733:   DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity

735:   Collective

737:   Input Parameter:
738: + dm - The DMMoab object
739: . vhandle - Vertex entity handle
740: . nconn - Number of entities whose coordinates are needed
741: - conn - The vertex entity handles

743:   Level: beginner

745: .seealso: DMMoabGetVertexCoordinates(), DMMoabGetVertexConnectivity()
746: @*/
747: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt* nconn, moab::EntityHandle **conn)
748: {
749:   PetscErrorCode  ierr;


755:   if (conn) {
756:     PetscFree(*conn);
757:   }
758:   if (nconn) *nconn = 0;
759:   return(0);
760: }


763: /*@
764:   DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity

766:   Collective

768:   Input Parameter:
769: + dm - The DMMoab object
770: - ehandle - Vertex entity handle

772:   Output Parameter:
773: + nconn - Number of entities whose coordinates are needed
774: - conn - The vertex entity handles

776:   Level: beginner

778: .seealso: DMMoabGetVertexCoordinates(), DMMoabGetVertexConnectivity(), DMMoabRestoreVertexConnectivity()
779: @*/
780: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt* nconn, const moab::EntityHandle **conn)
781: {
782:   DM_Moab        *dmmoab;
783:   const moab::EntityHandle *connect;
784:   std::vector<moab::EntityHandle> vconn;
785:   moab::ErrorCode merr;
786:   PetscInt nnodes;

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

793:   /* Get connectivity information in MOAB canonical ordering */
794:   merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes); MBERRNM(merr);
795:   if (conn) *conn = connect;
796:   if (nconn) *nconn = nnodes;
797:   return(0);
798: }


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

804:   Collective

806:   Input Parameter:
807: + dm - The DMMoab object
808: - ent - Entity handle

810:   Output Parameter:
811: . ent_on_boundary - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise

813:   Level: beginner

815: .seealso: DMMoabCheckBoundaryVertices()
816: @*/
817: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool* ent_on_boundary)
818: {
819:   moab::EntityType etype;
820:   DM_Moab         *dmmoab;
821:   PetscInt         edim;

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

828:   /* get the entity type and handle accordingly */
829:   etype = dmmoab->mbiface->type_from_handle(ent);
830:   if (etype >= moab::MBPOLYHEDRON) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Entity type on the boundary skin is invalid. EntityType = %D\n", etype);

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

835:   *ent_on_boundary = PETSC_FALSE;
836:   if (etype == moab::MBVERTEX && edim == 0) {
837:     *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
838:   }
839:   else {
840:     if (edim == dmmoab->dim) {  /* check the higher-dimensional elements first */
841:       if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
842:     }
843:     else {                      /* next check the lower-dimensional faces */
844:       if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
845:     }
846:   }
847:   return(0);
848: }


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

854:   Input Parameter:
855: + dm - The DMMoab object
856: . nconn - Number of handles
857: - cnt - Array of entity handles

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

862:   Level: beginner

864: .seealso: DMMoabIsEntityOnBoundary()
865: @*/
866: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool* isbdvtx)
867: {
868:   DM_Moab        *dmmoab;
869:   PetscInt       i;

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

877:   for (i = 0; i < nconn; ++i) {
878:     isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
879:   }
880:   return(0);
881: }


884: /*@
885:   DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary

887:   Input Parameter:
888: . dm - The DMMoab object

890:   Output Parameter:
891: + bdvtx - Boundary vertices
892: . bdelems - Boundary elements
893: - bdfaces - Boundary faces

895:   Level: beginner

897: .seealso: DMMoabCheckBoundaryVertices(), DMMoabIsEntityOnBoundary()
898: @*/
899: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range** bdelems, const moab::Range** bdfaces)
900: {
901:   DM_Moab        *dmmoab;

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

907:   if (bdvtx)  *bdvtx = dmmoab->bndyvtx;
908:   if (bdfaces)  *bdfaces = dmmoab->bndyfaces;
909:   if (bdelems)  *bdfaces = dmmoab->bndyelems;
910:   return(0);
911: }


914: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
915: {
916:   PetscErrorCode  ierr;
917:   PetscInt        i;
918:   moab::ErrorCode merr;
919:   DM_Moab        *dmmoab = (DM_Moab*)dm->data;


924:   dmmoab->refct--;
925:   if (!dmmoab->refct) {
926:     delete dmmoab->vlocal;
927:     delete dmmoab->vowned;
928:     delete dmmoab->vghost;
929:     delete dmmoab->elocal;
930:     delete dmmoab->eghost;
931:     delete dmmoab->bndyvtx;
932:     delete dmmoab->bndyfaces;
933:     delete dmmoab->bndyelems;

935:     PetscFree(dmmoab->gsindices);
936:     PetscFree2(dmmoab->gidmap, dmmoab->lidmap);
937:     PetscFree(dmmoab->dfill);
938:     PetscFree(dmmoab->ofill);
939:     PetscFree(dmmoab->materials);
940:     if (dmmoab->fieldNames) {
941:       for (i = 0; i < dmmoab->numFields; i++) {
942:         PetscFree(dmmoab->fieldNames[i]);
943:       }
944:       PetscFree(dmmoab->fieldNames);
945:     }

947:     if (dmmoab->nhlevels) {
948:       PetscFree(dmmoab->hsets);
949:       dmmoab->nhlevels = 0;
950:       if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
951:       dmmoab->hierarchy = NULL;
952:     }

954:     if (dmmoab->icreatedinstance) {
955:       delete dmmoab->pcomm;
956:       merr = dmmoab->mbiface->delete_mesh(); MBERRNM(merr);
957:       delete dmmoab->mbiface;
958:     }
959:     dmmoab->mbiface = NULL;
960: #ifdef MOAB_HAVE_MPI
961:     dmmoab->pcomm = NULL;
962: #endif
963:     VecScatterDestroy(&dmmoab->ltog_sendrecv);
964:     ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map);
965:     PetscFree(dm->data);
966:   }
967:   return(0);
968: }


971: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(PetscOptionItems *PetscOptionsObject, DM dm)
972: {
974:   DM_Moab        *dmmoab = (DM_Moab*)dm->data;

978:   PetscOptionsHead(PetscOptionsObject, "DMMoab Options");
979:   PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL,0);
980:   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);
981:   /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
982:   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, PETSC_MAX_PATH_LEN, NULL);
983:   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, PETSC_MAX_PATH_LEN, NULL);
984:   PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum*)&dmmoab->read_mode, NULL);
985:   PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum*)&dmmoab->write_mode, NULL);
986:   return(0);
987: }


990: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
991: {
992:   PetscErrorCode          ierr;
993:   moab::ErrorCode         merr;
994:   Vec                     local, global;
995:   IS                      from, to;
996:   moab::Range::iterator   iter;
997:   PetscInt                i, j, f, bs, vent, totsize, *lgmap;
998:   DM_Moab                *dmmoab = (DM_Moab*)dm->data;
999:   moab::Range             adjs;

1003:   /* Get the local and shared vertices and cache it */
1004:   if (dmmoab->mbiface == NULL) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB Interface before calling SetUp.");
1005: #ifdef MOAB_HAVE_MPI
1006:   if (dmmoab->pcomm == NULL) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB ParallelComm object before calling SetUp.");
1007: #endif

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

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

1019:     /* filter all the non-owned and shared entities out of the list */
1020:     // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
1021:     adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
1022:     merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost); MBERRNM(merr);
1023:     adjs = moab::subtract(adjs, *dmmoab->vghost);
1024:     *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
1025: #else
1026:     *dmmoab->vowned = *dmmoab->vlocal;
1027: #endif

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

1033: #ifdef MOAB_HAVE_MPI
1034:     MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
1035:     PetscInfo4(NULL, "Filset ID: %u, Vertices: local - %D, owned - %D, ghosted - %D.\n", dmmoab->fileset, dmmoab->vlocal->size(), dmmoab->nloc, dmmoab->nghost);
1036: #else
1037:     dmmoab->n = dmmoab->nloc;
1038: #endif
1039:   }

1041:   {
1042:     /* get the information about the local elements in the mesh */
1043:     dmmoab->eghost->clear();

1045:     /* first decipher the leading dimension */
1046:     for (i = 3; i > 0; i--) {
1047:       dmmoab->elocal->clear();
1048:       merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false); MBERRNM(merr);

1050:       /* store the current mesh dimension */
1051:       if (dmmoab->elocal->size()) {
1052:         dmmoab->dim = i;
1053:         break;
1054:       }
1055:     }

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

1059: #ifdef MOAB_HAVE_MPI
1060:     /* filter the ghosted and owned element list */
1061:     *dmmoab->eghost = *dmmoab->elocal;
1062:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1063:     *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
1064: #endif

1066:     dmmoab->neleloc = dmmoab->elocal->size();
1067:     dmmoab->neleghost = dmmoab->eghost->size();

1069: #ifdef MOAB_HAVE_MPI
1070:     MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
1071:     PetscInfo3(NULL, "%d-dim elements: owned - %D, ghosted - %D.\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost);
1072: #else
1073:     dmmoab->nele = dmmoab->neleloc;
1074: #endif
1075:   }

1077:   bs = dmmoab->bs;
1078:   if (!dmmoab->ltog_tag) {
1079:     /* Get the global ID tag. The global ID tag is applied to each
1080:        vertex. It acts as an global identifier which MOAB uses to
1081:        assemble the individual pieces of the mesh */
1082:     merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag); MBERRNM(merr);
1083:   }

1085:   totsize = dmmoab->vlocal->size();
1086:   if (totsize != dmmoab->nloc + dmmoab->nghost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Mismatch between local and owned+ghost vertices. %D != %D.", totsize, dmmoab->nloc + dmmoab->nghost);
1087:   PetscCalloc1(totsize, &dmmoab->gsindices);
1088:   {
1089:     /* first get the local indices */
1090:     merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]); MBERRNM(merr);
1091:     if (dmmoab->nghost) {  /* next get the ghosted indices */
1092:       merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]); MBERRNM(merr);
1093:     }

1095:     /* find out the local and global minima of GLOBAL_ID */
1096:     dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1097:     for (i = 0; i < totsize; ++i) {
1098:       if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1099:       if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1100:     }

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

1105:     /* set the GID map */
1106:     for (i = 0; i < totsize; ++i) {
1107:       dmmoab->gsindices[i] -= dmmoab->gminmax[0]; /* zero based index needed for IS */

1109:     }
1110:     dmmoab->lminmax[0] -= dmmoab->gminmax[0];
1111:     dmmoab->lminmax[1] -= dmmoab->gminmax[0];

1113:     PetscInfo4(NULL, "GLOBAL_ID: Local [min, max] - [%D, %D], Global [min, max] - [%D, %D]\n", dmmoab->lminmax[0], dmmoab->lminmax[1], dmmoab->gminmax[0], dmmoab->gminmax[1]);
1114:   }
1115:   if (!(dmmoab->bs == dmmoab->numFields || dmmoab->bs == 1)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mismatch between block size and number of component fields. %D != 1 OR %D != %D.", dmmoab->bs, dmmoab->bs, dmmoab->numFields);

1117:   {
1118:     dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1119:     dmmoab->seqend = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1120:     PetscInfo2(NULL, "SEQUENCE: Local [min, max] - [%D, %D]\n", dmmoab->seqstart, dmmoab->seqend);

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

1125:     i = j = 0;
1126:     /* set the owned vertex data first */
1127:     for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1128:       vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1129:       dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1130:       dmmoab->lidmap[vent] = i;
1131:       for (f = 0; f < dmmoab->numFields; f++, j++) {
1132:         lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1133:       }
1134:     }
1135:     /* next arrange all the ghosted data information */
1136:     for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1137:       vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1138:       dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1139:       dmmoab->lidmap[vent] = i;
1140:       for (f = 0; f < dmmoab->numFields; f++, j++) {
1141:         lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1142:       }
1143:     }

1145:     /* We need to create the Global to Local Vector Scatter Contexts
1146:        1) First create a local and global vector
1147:        2) Create a local and global IS
1148:        3) Create VecScatter and LtoGMapping objects
1149:        4) Cleanup the IS and Vec objects
1150:     */
1151:     DMCreateGlobalVector(dm, &global);
1152:     DMCreateLocalVector(dm, &local);

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

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

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

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

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

1172:     /* clean up IS, Vec */
1173:     PetscFree(lgmap);
1174:     ISDestroy(&from);
1175:     ISDestroy(&to);
1176:     VecDestroy(&local);
1177:     VecDestroy(&global);
1178:   }

1180:   dmmoab->bndyvtx = new moab::Range();
1181:   dmmoab->bndyfaces = new moab::Range();
1182:   dmmoab->bndyelems = new moab::Range();
1183:   /* skin the boundary and store nodes */
1184:   if (!dmmoab->hlevel) {
1185:     /* get the skin vertices of boundary faces for the current partition and then filter
1186:        the local, boundary faces, vertices and elements alone via PSTATUS flags;
1187:        this should not give us any ghosted boundary, but if user needs such a functionality
1188:        it would be easy to add it based on the find_skin query below */
1189:     moab::Skinner skinner(dmmoab->mbiface);

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

1194: #ifdef MOAB_HAVE_MPI
1195:     /* filter all the non-owned and shared entities out of the list */
1196:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1197:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT); MBERRNM(merr);
1198: #endif

1200:     /* get all the nodes via connectivity and the parent elements via adjacency information */
1201:     merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false); MBERRNM(ierr);
1202:     merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION); MBERRNM(ierr);
1203:   }
1204:   else {
1205:     /* Let us query the hierarchy manager and get the results directly for this level */
1206:     for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1207:       moab::EntityHandle elemHandle = *iter;
1208:       if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1209:         dmmoab->bndyelems->insert(elemHandle);
1210:         /* For this boundary element, query the vertices and add them to the list */
1211:         std::vector<moab::EntityHandle> connect;
1212:         merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect); MBERRNM(ierr);
1213:         for (unsigned iv=0; iv < connect.size(); ++iv)
1214:           if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv]))
1215:             dmmoab->bndyvtx->insert(connect[iv]);
1216:         /* Next, let us query the boundary faces and add them also to the list */
1217:         std::vector<moab::EntityHandle> faces;
1218:         merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim-1, faces); MBERRNM(ierr);
1219:         for (unsigned ifa=0; ifa < faces.size(); ++ifa)
1220:           if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa]))
1221:             dmmoab->bndyfaces->insert(faces[ifa]);
1222:       }
1223:     }
1224: #ifdef MOAB_HAVE_MPI
1225:     /* filter all the non-owned and shared entities out of the list */
1226:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx,   PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1227:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1228:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1229: #endif

1231:   }
1232:   PetscInfo3(NULL, "Found %D boundary vertices, %D boundary faces and %D boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size());

1234:   /* Get the material sets and populate the data for all locally owned elements */
1235:   {
1236:     PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials);
1237:     /* Get the count of entities of particular type from dmmoab->elocal
1238:        -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1239:     moab::Range msets;
1240:     merr = dmmoab->mbiface->get_entities_by_type_and_tag(dmmoab->fileset, moab::MBENTITYSET, &dmmoab->material_tag, NULL, 1, msets, moab::Interface::UNION);MB_CHK_ERR(merr);
1241:     if (msets.size() == 0) {
1242:       PetscInfo(NULL, "No material sets found in the fileset.");
1243:     }

1245:     for (unsigned i=0; i < msets.size(); ++i) {
1246:       moab::Range msetelems;
1247:       merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);MB_CHK_ERR(merr);
1248: #ifdef MOAB_HAVE_MPI
1249:       /* filter all the non-owned and shared entities out of the list */
1250:       merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1251: #endif

1253:       int partID;
1254:       moab::EntityHandle mset=msets[i];
1255:       merr = dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID);MB_CHK_ERR(merr);

1257:       for (unsigned j=0; j < msetelems.size(); ++j)
1258:         dmmoab->materials[dmmoab->elocal->index(msetelems[j])]=partID;
1259:     }
1260:   }

1262:   return(0);
1263: }


1266: /*@
1267:   DMMoabCreateVertices - Creates and adds several vertices to the primary set represented by the DM.

1269:   Collective

1271:   Input Parameters:
1272: + dm - The DM object
1273: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1274: . conn - The connectivity of the element
1275: - nverts - The number of vertices that form the element

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

1280:   Level: beginner

1282: .seealso: DMMoabCreateSubmesh(), DMMoabCreateElement()
1283: @*/
1284: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal* coords, PetscInt nverts, moab::Range* overts)
1285: {
1286:   moab::ErrorCode     merr;
1287:   DM_Moab            *dmmoab;
1288:   moab::Range         verts;


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

1296:   /* Insert new points */
1297:   merr = dmmoab->mbiface->create_vertices(&coords[0], nverts, verts); MBERRNM(merr);
1298:   merr = dmmoab->mbiface->add_entities(dmmoab->fileset, verts); MBERRNM(merr);

1300:   if (overts) *overts = verts;
1301:   return(0);
1302: }


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

1308:   Collective

1310:   Input Parameters:
1311: + dm - The DM object
1312: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1313: . conn - The connectivity of the element
1314: - nverts - The number of vertices that form the element

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

1319:   Level: beginner

1321: .seealso: DMMoabCreateSubmesh(), DMMoabCreateVertices()
1322: @*/
1323: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle* conn, PetscInt nverts, moab::EntityHandle* oelem)
1324: {
1325:   moab::ErrorCode     merr;
1326:   DM_Moab            *dmmoab;
1327:   moab::EntityHandle  elem;


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

1335:   /* Insert new element */
1336:   merr = dmmoab->mbiface->create_element(type, conn, nverts, elem); MBERRNM(merr);
1337:   merr = dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1); MBERRNM(merr);

1339:   if (oelem) *oelem = elem;
1340:   return(0);
1341: }


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

1349:   Collective

1351:   Input Parameters:
1352: . dm - The DM object

1354:   Output Parameter:
1355: . newdm  - The sub DM object with updated set information

1357:   Level: advanced

1359: .seealso: DMCreate(), DMMoabCreateVertices(), DMMoabCreateElement()
1360: @*/
1361: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1362: {
1363:   DM_Moab            *dmmoab;
1364:   DM_Moab            *ndmmoab;
1365:   moab::ErrorCode    merr;
1366:   PetscErrorCode     ierr;


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

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

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

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

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

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

1389:   /* preserve the field association between the parent and sub-mesh objects */
1390:   DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames);
1391:   return(0);
1392: }


1395: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1396: {
1397:   DM_Moab          *dmmoab = (DM_Moab*)(dm)->data;
1398:   const char       *name;
1399:   MPI_Comm          comm;
1400:   PetscMPIInt       size;
1401:   PetscErrorCode    ierr;

1404:   PetscObjectGetComm((PetscObject)dm, &comm);
1405:   MPI_Comm_size(comm, &size);
1406:   PetscObjectGetName((PetscObject) dm, &name);
1407:   PetscViewerASCIIPushTab(viewer);
1408:   if (name) {PetscViewerASCIIPrintf(viewer, "%s in %D dimensions:\n", name, dmmoab->dim);}
1409:   else      {PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dmmoab->dim);}
1410:   /* print details about the global mesh */
1411:   {
1412:     PetscViewerASCIIPushTab(viewer);
1413:     PetscViewerASCIIPrintf(viewer, "Sizes: cells=%D, vertices=%D, blocks=%D\n", dmmoab->nele, dmmoab->n, dmmoab->bs);
1414:     /* print boundary data */
1415:     PetscViewerASCIIPrintf(viewer, "Boundary trace:\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size());
1416:     {
1417:       PetscViewerASCIIPushTab(viewer);
1418:       PetscViewerASCIIPrintf(viewer, "cells=%D, faces=%D, vertices=%D\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size());
1419:       PetscViewerASCIIPopTab(viewer);
1420:     }
1421:     /* print field data */
1422:     PetscViewerASCIIPrintf(viewer, "Fields: %D components\n", dmmoab->numFields);
1423:     {
1424:       PetscViewerASCIIPushTab(viewer);
1425:       for (int i = 0; i < dmmoab->numFields; ++i) {
1426:         PetscViewerASCIIPrintf(viewer, "[%D] - %s\n", i, dmmoab->fieldNames[i]);
1427:       }
1428:       PetscViewerASCIIPopTab(viewer);
1429:     }
1430:     PetscViewerASCIIPopTab(viewer);
1431:   }
1432:   PetscViewerASCIIPopTab(viewer);
1433:   PetscViewerFlush(viewer);
1434:   return(0);
1435: }

1437: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1438: {
1439:   return(0);
1440: }

1442: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1443: {
1444:   return(0);
1445: }

1447: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1448: {
1449:   PetscBool      iascii, ishdf5, isvtk;

1455:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1456:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK,   &isvtk);
1457:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5,  &ishdf5);
1458:   if (iascii) {
1459:     DMMoabView_Ascii(dm, viewer);
1460:   } else if (ishdf5) {
1461: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1462:     PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
1463:     DMMoabView_HDF5(dm, viewer);
1464:     PetscViewerPopFormat(viewer);
1465: #else
1466:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1467: #endif
1468:   }
1469:   else if (isvtk) {
1470:     DMMoabView_VTK(dm, viewer);
1471:   }
1472:   return(0);
1473: }


1476: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
1477: {
1479:   dm->ops->view                            = DMView_Moab;
1480:   dm->ops->load                            = NULL /* DMLoad_Moab */;
1481:   dm->ops->setfromoptions                  = DMSetFromOptions_Moab;
1482:   dm->ops->clone                           = DMClone_Moab;
1483:   dm->ops->setup                           = DMSetUp_Moab;
1484:   dm->ops->createdefaultsection            = NULL;
1485:   dm->ops->createdefaultconstraints        = NULL;
1486:   dm->ops->createglobalvector              = DMCreateGlobalVector_Moab;
1487:   dm->ops->createlocalvector               = DMCreateLocalVector_Moab;
1488:   dm->ops->getlocaltoglobalmapping         = NULL;
1489:   dm->ops->createfieldis                   = NULL;
1490:   dm->ops->createcoordinatedm              = NULL /* DMCreateCoordinateDM_Moab */;
1491:   dm->ops->getcoloring                     = NULL;
1492:   dm->ops->creatematrix                    = DMCreateMatrix_Moab;
1493:   dm->ops->createinterpolation             = DMCreateInterpolation_Moab;
1494:   dm->ops->createinjection                 = NULL /* DMCreateInjection_Moab */;
1495:   dm->ops->refine                          = DMRefine_Moab;
1496:   dm->ops->coarsen                         = DMCoarsen_Moab;
1497:   dm->ops->refinehierarchy                 = DMRefineHierarchy_Moab;
1498:   dm->ops->coarsenhierarchy                = DMCoarsenHierarchy_Moab;
1499:   dm->ops->globaltolocalbegin              = DMGlobalToLocalBegin_Moab;
1500:   dm->ops->globaltolocalend                = DMGlobalToLocalEnd_Moab;
1501:   dm->ops->localtoglobalbegin              = DMLocalToGlobalBegin_Moab;
1502:   dm->ops->localtoglobalend                = DMLocalToGlobalEnd_Moab;
1503:   dm->ops->destroy                         = DMDestroy_Moab;
1504:   dm->ops->createsubdm                     = NULL /* DMCreateSubDM_Moab */;
1505:   dm->ops->getdimpoints                    = NULL /* DMGetDimPoints_Moab */;
1506:   dm->ops->locatepoints                    = NULL /* DMLocatePoints_Moab */;
1507:   return(0);
1508: }


1511: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1512: {
1513:   PetscErrorCode     ierr;

1516:   /* get all the necessary handles from the private DM object */
1517:   (*newdm)->data = (DM_Moab*) dm->data;
1518:   ((DM_Moab*)dm->data)->refct++;

1520:   PetscObjectChangeTypeName((PetscObject) *newdm, DMMOAB);
1521:   DMInitialize_Moab(*newdm);
1522:   return(0);
1523: }


1526: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1527: {

1532:   PetscNewLog(dm, (DM_Moab**)&dm->data);

1534:   ((DM_Moab*)dm->data)->bs = 1;
1535:   ((DM_Moab*)dm->data)->numFields = 1;
1536:   ((DM_Moab*)dm->data)->n = 0;
1537:   ((DM_Moab*)dm->data)->nloc = 0;
1538:   ((DM_Moab*)dm->data)->nghost = 0;
1539:   ((DM_Moab*)dm->data)->nele = 0;
1540:   ((DM_Moab*)dm->data)->neleloc = 0;
1541:   ((DM_Moab*)dm->data)->neleghost = 0;
1542:   ((DM_Moab*)dm->data)->ltog_map = NULL;
1543:   ((DM_Moab*)dm->data)->ltog_sendrecv = NULL;

1545:   ((DM_Moab*)dm->data)->refct = 1;
1546:   ((DM_Moab*)dm->data)->parent = NULL;
1547:   ((DM_Moab*)dm->data)->vlocal = new moab::Range();
1548:   ((DM_Moab*)dm->data)->vowned = new moab::Range();
1549:   ((DM_Moab*)dm->data)->vghost = new moab::Range();
1550:   ((DM_Moab*)dm->data)->elocal = new moab::Range();
1551:   ((DM_Moab*)dm->data)->eghost = new moab::Range();

1553:   DMInitialize_Moab(dm);
1554:   return(0);
1555: }