Actual source code: dmmoab.cxx

petsc-3.11.4 2019-09-28
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: http://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 on MPI_Comm

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

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

 64:   Level: beginner

 66: .keywords: DMMoab, create
 67: @*/
 68: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
 69: {

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

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

 82:   Collective on MPI_Comm

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

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

 95:   Level: intermediate

 97: .keywords: DMMoab, create
 98: @*/
 99: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
100: {
102:   moab::ErrorCode merr;
103:   DM             dmmb;
104:   DM_Moab        *dmmoab;


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

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

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

126: #ifdef MOAB_HAVE_MPI
127:   moab::EntityHandle partnset;

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

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

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

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

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

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


169: #ifdef MOAB_HAVE_MPI

171: /*@
172:   DMMoabGetParallelComm - Get the ParallelComm used with this DMMoab

174:   Collective on MPI_Comm

176:   Input Parameter:
177: . dm    - The DMMoab object being set

179:   Output Parameter:
180: . pcomm - The ParallelComm for the DMMoab

182:   Level: beginner

184: .keywords: DMMoab, create
185: @*/
186: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
187: {
190:   *pcomm = ((DM_Moab*)(dm)->data)->pcomm;
191:   return(0);
192: }

194: #endif /* MOAB_HAVE_MPI */


197: /*@
198:   DMMoabSetInterface - Set the MOAB instance used with this DMMoab

200:   Collective on MPI_Comm

202:   Input Parameter:
203: . dm      - The DMMoab object being set
204: . mbiface - The MOAB instance being set on this DMMoab

206:   Level: beginner

208: .keywords: DMMoab, create
209: @*/
210: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
211: {
212:   DM_Moab        *dmmoab = (DM_Moab*)(dm)->data;

217: #ifdef MOAB_HAVE_MPI
218:   dmmoab->pcomm = NULL;
219: #endif
220:   dmmoab->mbiface = mbiface;
221:   dmmoab->icreatedinstance = PETSC_FALSE;
222:   return(0);
223: }


226: /*@
227:   DMMoabGetInterface - Get the MOAB instance used with this DMMoab

229:   Collective on MPI_Comm

231:   Input Parameter:
232: . dm      - The DMMoab object being set

234:   Output Parameter:
235: . mbiface - The MOAB instance set on this DMMoab

237:   Level: beginner

239: .keywords: DMMoab, create
240: @*/
241: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
242: {
243:   PetscErrorCode   ierr;
244:   static PetscBool cite = PETSC_FALSE;

248:   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);
249:   *mbiface = ((DM_Moab*)dm->data)->mbiface;
250:   return(0);
251: }


254: /*@
255:   DMMoabSetLocalVertices - Set the entities having DOFs on this DMMoab

257:   Collective on MPI_Comm

259:   Input Parameter:
260: . dm    - The DMMoab object being set
261: . range - The entities treated by this DMMoab

263:   Level: beginner

265: .keywords: DMMoab, create
266: @*/
267: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
268: {
269:   moab::Range     tmpvtxs;
270:   DM_Moab        *dmmoab = (DM_Moab*)(dm)->data;

274:   dmmoab->vlocal->clear();
275:   dmmoab->vowned->clear();

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

279: #ifdef MOAB_HAVE_MPI
280:   moab::ErrorCode merr;
281:   /* filter based on parallel status */
282:   merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned); MBERRNM(merr);

284:   /* filter all the non-owned and shared entities out of the list */
285:   tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
286:   merr = dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost); MBERRNM(merr);
287:   tmpvtxs = moab::subtract(tmpvtxs, *dmmoab->vghost);
288:   *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
289: #else
290:   *dmmoab->vowned = *dmmoab->vlocal;
291: #endif

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


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

309:   Collective on MPI_Comm

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

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

317:   Level: beginner

319: .keywords: DMMoab, create
320: @*/
321: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
322: {
325:   if (local) *local = *((DM_Moab*)dm->data)->vlocal;
326:   return(0);
327: }



331: /*@
332:   DMMoabGetLocalVertices - Get the entities having DOFs on this DMMoab

334:   Collective on MPI_Comm

336:   Input Parameter:
337: . dm    - The DMMoab object being set

339:   Output Parameter:
340: . owned - The owned vertex entities in this DMMoab
341: . ghost - The ghosted entities (non-owned) stored locally in this partition

343:   Level: beginner

345: .keywords: DMMoab, create
346: @*/
347: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
348: {
351:   if (owned) *owned = ((DM_Moab*)dm->data)->vowned;
352:   if (ghost) *ghost = ((DM_Moab*)dm->data)->vghost;
353:   return(0);
354: }

356: /*@
357:   DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned

359:   Collective on MPI_Comm

361:   Input Parameter:
362: . dm    - The DMMoab object being set

364:   Output Parameter:
365: . range - The entities owned locally

367:   Level: beginner

369: .keywords: DMMoab, create
370: @*/
371: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
372: {
375:   if (range) *range = ((DM_Moab*)dm->data)->elocal;
376:   return(0);
377: }


380: /*@
381:   DMMoabSetLocalElements - Set the entities having DOFs on this DMMoab

383:   Collective on MPI_Comm

385:   Input Parameter:
386: . dm    - The DMMoab object being set
387: . range - The entities treated by this DMMoab

389:   Level: beginner

391: .keywords: DMMoab, create
392: @*/
393: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
394: {
395:   DM_Moab        *dmmoab = (DM_Moab*)(dm)->data;

399:   dmmoab->elocal->clear();
400:   dmmoab->eghost->clear();
401:   dmmoab->elocal->insert(range->begin(), range->end());
402: #ifdef MOAB_HAVE_MPI
403:   moab::ErrorCode merr;
404:   merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
405:   *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
406: #endif
407:   dmmoab->neleloc = dmmoab->elocal->size();
408:   dmmoab->neleghost = dmmoab->eghost->size();
409: #ifdef MOAB_HAVE_MPI
410:   PetscErrorCode  ierr;
411:   MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
412:   PetscInfo2(dm, "Created %D local and %D global elements.\n", dmmoab->neleloc, dmmoab->nele);
413: #else
414:   dmmoab->nele = dmmoab->neleloc;
415: #endif
416:   return(0);
417: }


420: /*@
421:   DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering

423:   Collective on MPI_Comm

425:   Input Parameter:
426: . dm      - The DMMoab object being set
427: . ltogtag - The MOAB tag used for local to global ids

429:   Level: beginner

431: .keywords: DMMoab, create
432: @*/
433: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
434: {
437:   ((DM_Moab*)dm->data)->ltog_tag = ltogtag;
438:   return(0);
439: }


442: /*@
443:   DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering

445:   Collective on MPI_Comm

447:   Input Parameter:
448: . dm      - The DMMoab object being set

450:   Output Parameter:
451: . ltogtag - The MOAB tag used for local to global ids

453:   Level: beginner

455: .keywords: DMMoab, create
456: @*/
457: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
458: {
461:   *ltog_tag = ((DM_Moab*)dm->data)->ltog_tag;
462:   return(0);
463: }


466: /*@
467:   DMMoabSetBlockSize - Set the block size used with this DMMoab

469:   Collective on MPI_Comm

471:   Input Parameter:
472: . dm - The DMMoab object being set
473: . bs - The block size used with this DMMoab

475:   Level: beginner

477: .keywords: DMMoab, create
478: @*/
479: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
480: {
483:   ((DM_Moab*)dm->data)->bs = bs;
484:   return(0);
485: }


488: /*@
489:   DMMoabGetBlockSize - Get the block size used with this DMMoab

491:   Collective on MPI_Comm

493:   Input Parameter:
494: . dm - The DMMoab object being set

496:   Output Parameter:
497: . bs - The block size used with this DMMoab

499:   Level: beginner

501: .keywords: DMMoab, create
502: @*/
503: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
504: {
507:   *bs = ((DM_Moab*)dm->data)->bs;
508:   return(0);
509: }


512: /*@
513:   DMMoabGetSize - Get the global vertex size used with this DMMoab

515:   Collective on DM

517:   Input Parameter:
518: . dm - The DMMoab object being set

520:   Output Parameter:
521: . neg - The number of global elements in the DMMoab instance
522: . nvg - The number of global vertices in the DMMoab instance

524:   Level: beginner

526: .keywords: DMMoab, create
527: @*/
528: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
529: {
532:   if (neg) *neg = ((DM_Moab*)dm->data)->nele;
533:   if (nvg) *nvg = ((DM_Moab*)dm->data)->n;
534:   return(0);
535: }


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

541:   Collective on DM

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

546:   Output Parameter:
547: + nel - The number of owned elements in this processor
548: . neg - The number of ghosted elements in this processor
549: . nvl - The number of owned vertices in this processor
550: . nvg - The number of ghosted vertices in this processor

552:   Level: beginner

554: .keywords: DMMoab, create
555: @*/
556: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
557: {
560:   if (nel) *nel = ((DM_Moab*)dm->data)->neleloc;
561:   if (neg) *neg = ((DM_Moab*)dm->data)->neleghost;
562:   if (nvl) *nvl = ((DM_Moab*)dm->data)->nloc;
563:   if (nvg) *nvg = ((DM_Moab*)dm->data)->nghost;
564:   return(0);
565: }


568: /*@
569:   DMMoabGetOffset - Get the local offset for the global vector

571:   Collective on MPI_Comm

573:   Input Parameter:
574: . dm - The DMMoab object being set

576:   Output Parameter:
577: . offset - The local offset for the global vector

579:   Level: beginner

581: .keywords: DMMoab, create
582: @*/
583: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
584: {
587:   *offset = ((DM_Moab*)dm->data)->vstart;
588:   return(0);
589: }


592: /*@
593:   DMMoabGetDimension - Get the dimension of the DM Mesh

595:   Collective on MPI_Comm

597:   Input Parameter:
598: . dm - The DMMoab object

600:   Output Parameter:
601: . dim - The dimension of DM

603:   Level: beginner

605: .keywords: DMMoab, create
606: @*/
607: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
608: {
611:   *dim = ((DM_Moab*)dm->data)->dim;
612:   return(0);
613: }


616: /*@
617:   DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
618:   generated through uniform refinement.

620:   Collective on DM

622:   Input Parameter:
623: . dm - The DMMoab object being set

625:   Output Parameter:
626: . nvg - The current mesh hierarchy level

628:   Level: beginner

630: .keywords: DMMoab, multigrid
631: @*/
632: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
633: {
636:   if (nlevel) *nlevel = ((DM_Moab*)dm->data)->hlevel;
637:   return(0);
638: }


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

644:   Collective on MPI_Comm

646:   Input Parameter:
647: . dm - The DMMoab object
648: . ehandle - The element entity handle

650:   Output Parameter:
651: . mat - The material ID for the current entity

653:   Level: beginner

655: .keywords: DMMoab, create
656: @*/
657: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
658: {
659:   DM_Moab         *dmmoab;

663:   if (*mat) {
664:     dmmoab = (DM_Moab*)(dm)->data;
665:     *mat = dmmoab->materials[dmmoab->elocal->index(ehandle)];
666:   }
667:   return(0);
668: }


671: /*@
672:   DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities

674:   Collective on MPI_Comm

676:   Input Parameter:
677: . dm - The DMMoab object
678: . nconn - Number of entities whose coordinates are needed
679: . conn - The vertex entity handles

681:   Output Parameter:
682: . vpos - The coordinates of the requested vertex entities

684:   Level: beginner

686: .seealso: DMMoabGetVertexConnectivity()
687: @*/
688: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
689: {
690:   DM_Moab         *dmmoab;
691:   moab::ErrorCode merr;

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

699:   /* Get connectivity information in MOAB canonical ordering */
700:   if (dmmoab->hlevel) {
701:     merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle*>(conn), nconn, dmmoab->hlevel, vpos);MBERRNM(merr);
702:   }
703:   else {
704:     merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);MBERRNM(merr);
705:   }
706:   return(0);
707: }


710: /*@
711:   DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity

713:   Collective on MPI_Comm

715:   Input Parameter:
716: . dm - The DMMoab object
717: . vhandle - Vertex entity handle

719:   Output Parameter:
720: . nconn - Number of entities whose coordinates are needed
721: . conn - The vertex entity handles

723:   Level: beginner

725: .seealso: DMMoabGetVertexCoordinates(), DMMoabRestoreVertexConnectivity()
726: @*/
727: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt* nconn, moab::EntityHandle **conn)
728: {
729:   DM_Moab        *dmmoab;
730:   std::vector<moab::EntityHandle> adj_entities, connect;
731:   PetscErrorCode  ierr;
732:   moab::ErrorCode merr;

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

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

743:   if (conn) {
744:     PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn);
745:     PetscMemcpy(*conn, &connect[0], sizeof(moab::EntityHandle) * connect.size());
746:   }
747:   if (nconn) *nconn = connect.size();
748:   return(0);
749: }


752: /*@
753:   DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity

755:   Collective on MPI_Comm

757:   Input Parameter:
758: . dm - The DMMoab object
759: . vhandle - Vertex entity handle
760: . nconn - Number of entities whose coordinates are needed
761: . conn - The vertex entity handles

763:   Level: beginner

765: .seealso: DMMoabGetVertexCoordinates(), DMMoabGetVertexConnectivity()
766: @*/
767: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt* nconn, moab::EntityHandle **conn)
768: {
769:   PetscErrorCode  ierr;


775:   if (conn) {
776:     PetscFree(*conn);
777:   }
778:   if (nconn) *nconn = 0;
779:   return(0);
780: }


783: /*@
784:   DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity

786:   Collective on MPI_Comm

788:   Input Parameter:
789: . dm - The DMMoab object
790: . ehandle - Vertex entity handle

792:   Output Parameter:
793: . nconn - Number of entities whose coordinates are needed
794: . conn - The vertex entity handles

796:   Level: beginner

798: .seealso: DMMoabGetVertexCoordinates(), DMMoabGetVertexConnectivity(), DMMoabRestoreVertexConnectivity()
799: @*/
800: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt* nconn, const moab::EntityHandle **conn)
801: {
802:   DM_Moab        *dmmoab;
803:   const moab::EntityHandle *connect;
804:   std::vector<moab::EntityHandle> vconn;
805:   moab::ErrorCode merr;
806:   PetscInt nnodes;

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

813:   /* Get connectivity information in MOAB canonical ordering */
814:   merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes); MBERRNM(merr);
815:   if (conn) *conn = connect;
816:   if (nconn) *nconn = nnodes;
817:   return(0);
818: }


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

824:   Collective on MPI_Comm

826:   Input Parameter:
827: . dm - The DMMoab object
828: . ent - Entity handle

830:   Output Parameter:
831: . ent_on_boundary - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise

833:   Level: beginner

835: .seealso: DMMoabCheckBoundaryVertices()
836: @*/
837: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool* ent_on_boundary)
838: {
839:   moab::EntityType etype;
840:   DM_Moab         *dmmoab;
841:   PetscInt         edim;

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

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

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

855:   *ent_on_boundary = PETSC_FALSE;
856:   if (etype == moab::MBVERTEX && edim == 0) {
857:     *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
858:   }
859:   else {
860:     if (edim == dmmoab->dim) {  /* check the higher-dimensional elements first */
861:       if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
862:     }
863:     else {                      /* next check the lower-dimensional faces */
864:       if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
865:     }
866:   }
867:   return(0);
868: }


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

874:   Input Parameter:
875: . dm - The DMMoab object
876: . nconn - Number of handles
877: . cnt - Array of entity handles

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

882:   Level: beginner

884: .seealso: DMMoabIsEntityOnBoundary()
885: @*/
886: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool* isbdvtx)
887: {
888:   DM_Moab        *dmmoab;
889:   PetscInt       i;

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

897:   for (i = 0; i < nconn; ++i) {
898:     isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
899:   }
900:   return(0);
901: }


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

907:   Input Parameter:
908: . dm - The DMMoab object

910:   Output Parameter:
911: . bdvtx - Boundary vertices
912: . bdelems - Boundary elements
913: . bdfaces - Boundary faces

915:   Level: beginner

917: .seealso: DMMoabCheckBoundaryVertices(), DMMoabIsEntityOnBoundary()
918: @*/
919: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range** bdelems, const moab::Range** bdfaces)
920: {
921:   DM_Moab        *dmmoab;

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

927:   if (bdvtx)  *bdvtx = dmmoab->bndyvtx;
928:   if (bdfaces)  *bdfaces = dmmoab->bndyfaces;
929:   if (bdelems)  *bdfaces = dmmoab->bndyelems;
930:   return(0);
931: }


934: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
935: {
936:   PetscErrorCode  ierr;
937:   PetscInt        i;
938:   moab::ErrorCode merr;
939:   DM_Moab        *dmmoab = (DM_Moab*)dm->data;


944:   dmmoab->refct--;
945:   if (!dmmoab->refct) {
946:     delete dmmoab->vlocal;
947:     delete dmmoab->vowned;
948:     delete dmmoab->vghost;
949:     delete dmmoab->elocal;
950:     delete dmmoab->eghost;
951:     delete dmmoab->bndyvtx;
952:     delete dmmoab->bndyfaces;
953:     delete dmmoab->bndyelems;

955:     PetscFree(dmmoab->gsindices);
956:     PetscFree2(dmmoab->gidmap, dmmoab->lidmap);
957:     PetscFree(dmmoab->dfill);
958:     PetscFree(dmmoab->ofill);
959:     PetscFree(dmmoab->materials);
960:     if (dmmoab->fieldNames) {
961:       for (i = 0; i < dmmoab->numFields; i++) {
962:         PetscFree(dmmoab->fieldNames[i]);
963:       }
964:       PetscFree(dmmoab->fieldNames);
965:     }

967:     if (dmmoab->nhlevels) {
968:       PetscFree(dmmoab->hsets);
969:       dmmoab->nhlevels = 0;
970:       if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
971:       dmmoab->hierarchy = NULL;
972:     }

974:     if (dmmoab->icreatedinstance) {
975:       delete dmmoab->pcomm;
976:       merr = dmmoab->mbiface->delete_mesh(); MBERRNM(merr);
977:       delete dmmoab->mbiface;
978:     }
979:     dmmoab->mbiface = NULL;
980: #ifdef MOAB_HAVE_MPI
981:     dmmoab->pcomm = NULL;
982: #endif
983:     VecScatterDestroy(&dmmoab->ltog_sendrecv);
984:     ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map);
985:     PetscFree(dm->data);
986:   }
987:   return(0);
988: }


991: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(PetscOptionItems *PetscOptionsObject, DM dm)
992: {
994:   DM_Moab        *dmmoab = (DM_Moab*)dm->data;

998:   PetscOptionsHead(PetscOptionsObject, "DMMoab Options");
999:   PetscOptionsInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL);
1000:   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);
1001:   /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
1002:   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);
1003:   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);
1004:   PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum*)&dmmoab->read_mode, NULL);
1005:   PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum*)&dmmoab->write_mode, NULL);
1006:   return(0);
1007: }


1010: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
1011: {
1012:   PetscErrorCode          ierr;
1013:   moab::ErrorCode         merr;
1014:   Vec                     local, global;
1015:   IS                      from, to;
1016:   moab::Range::iterator   iter;
1017:   PetscInt                i, j, f, bs, vent, totsize, *lgmap;
1018:   DM_Moab                *dmmoab = (DM_Moab*)dm->data;
1019:   moab::Range             adjs;

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

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

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

1039:     /* filter all the non-owned and shared entities out of the list */
1040:     // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
1041:     adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
1042:     merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost); MBERRNM(merr);
1043:     adjs = moab::subtract(adjs, *dmmoab->vghost);
1044:     *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
1045: #else
1046:     *dmmoab->vowned = *dmmoab->vlocal;
1047: #endif

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

1053: #ifdef MOAB_HAVE_MPI
1054:     MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
1055:     PetscInfo4(NULL, "Filset ID: %u, Vertices: local - %D, owned - %D, ghosted - %D.\n", dmmoab->fileset, dmmoab->vlocal->size(), dmmoab->nloc, dmmoab->nghost);
1056: #else
1057:     dmmoab->n = dmmoab->nloc;
1058: #endif
1059:   }

1061:   {
1062:     /* get the information about the local elements in the mesh */
1063:     dmmoab->eghost->clear();

1065:     /* first decipher the leading dimension */
1066:     for (i = 3; i > 0; i--) {
1067:       dmmoab->elocal->clear();
1068:       merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false); MBERRNM(merr);

1070:       /* store the current mesh dimension */
1071:       if (dmmoab->elocal->size()) {
1072:         dmmoab->dim = i;
1073:         break;
1074:       }
1075:     }

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

1079: #ifdef MOAB_HAVE_MPI
1080:     /* filter the ghosted and owned element list */
1081:     *dmmoab->eghost = *dmmoab->elocal;
1082:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1083:     *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
1084: #endif

1086:     dmmoab->neleloc = dmmoab->elocal->size();
1087:     dmmoab->neleghost = dmmoab->eghost->size();

1089: #ifdef MOAB_HAVE_MPI
1090:     MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
1091:     PetscInfo3(NULL, "%d-dim elements: owned - %D, ghosted - %D.\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost);
1092: #else
1093:     dmmoab->nele = dmmoab->neleloc;
1094: #endif
1095:   }

1097:   bs = dmmoab->bs;
1098:   if (!dmmoab->ltog_tag) {
1099:     /* Get the global ID tag. The global ID tag is applied to each
1100:        vertex. It acts as an global identifier which MOAB uses to
1101:        assemble the individual pieces of the mesh */
1102:     merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag); MBERRNM(merr);
1103:   }

1105:   totsize = dmmoab->vlocal->size();
1106:   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);
1107:   PetscCalloc1(totsize, &dmmoab->gsindices);
1108:   {
1109:     /* first get the local indices */
1110:     merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]); MBERRNM(merr);
1111:     if (dmmoab->nghost) {  /* next get the ghosted indices */
1112:       merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]); MBERRNM(merr);
1113:     }

1115:     /* find out the local and global minima of GLOBAL_ID */
1116:     dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1117:     for (i = 0; i < totsize; ++i) {
1118:       if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1119:       if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1120:     }

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

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

1129:     }
1130:     dmmoab->lminmax[0] -= dmmoab->gminmax[0];
1131:     dmmoab->lminmax[1] -= dmmoab->gminmax[0];

1133:     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]);
1134:   }
1135:   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);

1137:   {
1138:     dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1139:     dmmoab->seqend = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1140:     PetscInfo2(NULL, "SEQUENCE: Local [min, max] - [%D, %D]\n", dmmoab->seqstart, dmmoab->seqend);

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

1145:     i = j = 0;
1146:     /* set the owned vertex data first */
1147:     for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1148:       vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1149:       dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1150:       dmmoab->lidmap[vent] = i;
1151:       for (f = 0; f < dmmoab->numFields; f++, j++) {
1152:         lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1153:       }
1154:     }
1155:     /* next arrange all the ghosted data information */
1156:     for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1157:       vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1158:       dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1159:       dmmoab->lidmap[vent] = i;
1160:       for (f = 0; f < dmmoab->numFields; f++, j++) {
1161:         lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1162:       }
1163:     }

1165:     /* We need to create the Global to Local Vector Scatter Contexts
1166:        1) First create a local and global vector
1167:        2) Create a local and global IS
1168:        3) Create VecScatter and LtoGMapping objects
1169:        4) Cleanup the IS and Vec objects
1170:     */
1171:     DMCreateGlobalVector(dm, &global);
1172:     DMCreateLocalVector(dm, &local);

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

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

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

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

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

1192:     /* clean up IS, Vec */
1193:     PetscFree(lgmap);
1194:     ISDestroy(&from);
1195:     ISDestroy(&to);
1196:     VecDestroy(&local);
1197:     VecDestroy(&global);
1198:   }

1200:   dmmoab->bndyvtx = new moab::Range();
1201:   dmmoab->bndyfaces = new moab::Range();
1202:   dmmoab->bndyelems = new moab::Range();
1203:   /* skin the boundary and store nodes */
1204:   if (!dmmoab->hlevel) {
1205:     /* get the skin vertices of boundary faces for the current partition and then filter
1206:        the local, boundary faces, vertices and elements alone via PSTATUS flags;
1207:        this should not give us any ghosted boundary, but if user needs such a functionality
1208:        it would be easy to add it based on the find_skin query below */
1209:     moab::Skinner skinner(dmmoab->mbiface);

1211:     /* get the entities on the skin - only the faces */
1212:     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

1214: #ifdef MOAB_HAVE_MPI
1215:     /* filter all the non-owned and shared entities out of the list */
1216:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1217:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT); MBERRNM(merr);
1218: #endif

1220:     /* get all the nodes via connectivity and the parent elements via adjacency information */
1221:     merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false); MBERRNM(ierr);
1222:     merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION); MBERRNM(ierr);
1223:   }
1224:   else {
1225:     /* Let us query the hierarchy manager and get the results directly for this level */
1226:     for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1227:       moab::EntityHandle elemHandle = *iter;
1228:       if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1229:         dmmoab->bndyelems->insert(elemHandle);
1230:         /* For this boundary element, query the vertices and add them to the list */
1231:         std::vector<moab::EntityHandle> connect;
1232:         merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect); MBERRNM(ierr);
1233:         for (unsigned iv=0; iv < connect.size(); ++iv)
1234:           if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv]))
1235:             dmmoab->bndyvtx->insert(connect[iv]);
1236:         /* Next, let us query the boundary faces and add them also to the list */
1237:         std::vector<moab::EntityHandle> faces;
1238:         merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim-1, faces); MBERRNM(ierr);
1239:         for (unsigned ifa=0; ifa < faces.size(); ++ifa)
1240:           if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa]))
1241:             dmmoab->bndyfaces->insert(faces[ifa]);
1242:       }
1243:     }
1244: #ifdef MOAB_HAVE_MPI
1245:     /* filter all the non-owned and shared entities out of the list */
1246:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx,   PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1247:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1248:     merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1249: #endif

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

1254:   /* Get the material sets and populate the data for all locally owned elements */
1255:   {
1256:     PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials);
1257:     /* Get the count of entities of particular type from dmmoab->elocal
1258:        -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1259:     moab::Range msets;
1260:     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);
1261:     if (msets.size() == 0) {
1262:       PetscInfo(NULL, "No material sets found in the fileset.");
1263:     }

1265:     for (unsigned i=0; i < msets.size(); ++i) {
1266:       moab::Range msetelems;
1267:       merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);MB_CHK_ERR(merr);
1268: #ifdef MOAB_HAVE_MPI
1269:       /* filter all the non-owned and shared entities out of the list */
1270:       merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1271: #endif

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

1277:       for (unsigned j=0; j < msetelems.size(); ++j)
1278:         dmmoab->materials[dmmoab->elocal->index(msetelems[j])]=partID;
1279:     }
1280:   }

1282:   return(0);
1283: }


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

1289:   Collective on MPI_Comm

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

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

1300:   Level: beginner

1302: .keywords: DM, create vertices

1304: .seealso: DMMoabCreateSubmesh(), DMMoabCreateElement()
1305: @*/
1306: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal* coords, PetscInt nverts, moab::Range* overts)
1307: {
1308:   moab::ErrorCode     merr;
1309:   DM_Moab            *dmmoab;
1310:   moab::Range         verts;


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

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

1322:   if (overts) *overts = verts;
1323:   return(0);
1324: }


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

1330:   Collective on MPI_Comm

1332:   Input Parameters:
1333: + dm - The DM object
1334: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1335: . conn - The connectivity of the element
1336: . nverts - The number of vertices that form the element

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

1341:   Level: beginner

1343: .keywords: DM, create element

1345: .seealso: DMMoabCreateSubmesh(), DMMoabCreateVertices()
1346: @*/
1347: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle* conn, PetscInt nverts, moab::EntityHandle* oelem)
1348: {
1349:   moab::ErrorCode     merr;
1350:   DM_Moab            *dmmoab;
1351:   moab::EntityHandle  elem;


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

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

1363:   if (oelem) *oelem = elem;
1364:   return(0);
1365: }


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

1373:   Collective on MPI_Comm

1375:   Input Parameters:
1376: + dm - The DM object

1378:   Output Parameter:
1379: . newdm  - The sub DM object with updated set information

1381:   Level: advanced

1383: .keywords: DM, sub-DM

1385: .seealso: DMCreate(), DMMoabCreateVertices(), DMMoabCreateElement()
1386: @*/
1387: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1388: {
1389:   DM_Moab            *dmmoab;
1390:   DM_Moab            *ndmmoab;
1391:   moab::ErrorCode    merr;
1392:   PetscErrorCode     ierr;


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

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

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

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

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

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

1415:   /* preserve the field association between the parent and sub-mesh objects */
1416:   DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames);
1417:   return(0);
1418: }


1421: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1422: {
1423:   DM_Moab          *dmmoab = (DM_Moab*)(dm)->data;
1424:   const char       *name;
1425:   MPI_Comm          comm;
1426:   PetscMPIInt       size;
1427:   PetscErrorCode    ierr;

1430:   PetscObjectGetComm((PetscObject)dm, &comm);
1431:   MPI_Comm_size(comm, &size);
1432:   PetscObjectGetName((PetscObject) dm, &name);
1433:   PetscViewerASCIIPushTab(viewer);
1434:   if (name) {PetscViewerASCIIPrintf(viewer, "%s in %D dimensions:\n", name, dmmoab->dim);}
1435:   else      {PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dmmoab->dim);}
1436:   /* print details about the global mesh */
1437:   {
1438:     PetscViewerASCIIPushTab(viewer);
1439:     PetscViewerASCIIPrintf(viewer, "Sizes: cells=%D, vertices=%D, blocks=%D\n", dmmoab->nele, dmmoab->n, dmmoab->bs);
1440:     /* print boundary data */
1441:     PetscViewerASCIIPrintf(viewer, "Boundary trace:\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size());
1442:     {
1443:       PetscViewerASCIIPushTab(viewer);
1444:       PetscViewerASCIIPrintf(viewer, "cells=%D, faces=%D, vertices=%D\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size());
1445:       PetscViewerASCIIPopTab(viewer);
1446:     }
1447:     /* print field data */
1448:     PetscViewerASCIIPrintf(viewer, "Fields: %D components\n", dmmoab->numFields);
1449:     {
1450:       PetscViewerASCIIPushTab(viewer);
1451:       for (int i = 0; i < dmmoab->numFields; ++i) {
1452:         PetscViewerASCIIPrintf(viewer, "[%D] - %s\n", i, dmmoab->fieldNames[i]);
1453:       }
1454:       PetscViewerASCIIPopTab(viewer);
1455:     }
1456:     PetscViewerASCIIPopTab(viewer);
1457:   }
1458:   PetscViewerASCIIPopTab(viewer);
1459:   PetscViewerFlush(viewer);
1460:   return(0);
1461: }

1463: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1464: {
1465:   return(0);
1466: }

1468: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1469: {
1470:   return(0);
1471: }

1473: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1474: {
1475:   PetscBool      iascii, ishdf5, isvtk;

1481:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1482:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK,   &isvtk);
1483:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5,  &ishdf5);
1484:   if (iascii) {
1485:     DMMoabView_Ascii(dm, viewer);
1486:   } else if (ishdf5) {
1487: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1488:     PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
1489:     DMMoabView_HDF5(dm, viewer);
1490:     PetscViewerPopFormat(viewer);
1491: #else
1492:     SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1493: #endif
1494:   }
1495:   else if (isvtk) {
1496:     DMMoabView_VTK(dm, viewer);
1497:   }
1498:   return(0);
1499: }


1502: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
1503: {
1505:   dm->ops->view                            = DMView_Moab;
1506:   dm->ops->load                            = NULL /* DMLoad_Moab */;
1507:   dm->ops->setfromoptions                  = DMSetFromOptions_Moab;
1508:   dm->ops->clone                           = DMClone_Moab;
1509:   dm->ops->setup                           = DMSetUp_Moab;
1510:   dm->ops->createdefaultsection            = NULL;
1511:   dm->ops->createdefaultconstraints        = NULL;
1512:   dm->ops->createglobalvector              = DMCreateGlobalVector_Moab;
1513:   dm->ops->createlocalvector               = DMCreateLocalVector_Moab;
1514:   dm->ops->getlocaltoglobalmapping         = NULL;
1515:   dm->ops->createfieldis                   = NULL;
1516:   dm->ops->createcoordinatedm              = NULL /* DMCreateCoordinateDM_Moab */;
1517:   dm->ops->getcoloring                     = NULL;
1518:   dm->ops->creatematrix                    = DMCreateMatrix_Moab;
1519:   dm->ops->createinterpolation             = DMCreateInterpolation_Moab;
1520:   dm->ops->getaggregates                   = NULL;
1521:   dm->ops->getinjection                    = NULL /* DMCreateInjection_Moab */;
1522:   dm->ops->refine                          = DMRefine_Moab;
1523:   dm->ops->coarsen                         = DMCoarsen_Moab;
1524:   dm->ops->refinehierarchy                 = DMRefineHierarchy_Moab;
1525:   dm->ops->coarsenhierarchy                = DMCoarsenHierarchy_Moab;
1526:   dm->ops->globaltolocalbegin              = DMGlobalToLocalBegin_Moab;
1527:   dm->ops->globaltolocalend                = DMGlobalToLocalEnd_Moab;
1528:   dm->ops->localtoglobalbegin              = DMLocalToGlobalBegin_Moab;
1529:   dm->ops->localtoglobalend                = DMLocalToGlobalEnd_Moab;
1530:   dm->ops->destroy                         = DMDestroy_Moab;
1531:   dm->ops->createsubdm                     = NULL /* DMCreateSubDM_Moab */;
1532:   dm->ops->getdimpoints                    = NULL /* DMGetDimPoints_Moab */;
1533:   dm->ops->locatepoints                    = NULL /* DMLocatePoints_Moab */;
1534:   return(0);
1535: }


1538: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1539: {
1540:   PetscErrorCode     ierr;

1543:   PetscObjectChangeTypeName((PetscObject) * newdm, DMMOAB);

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

1549:   DMInitialize_Moab(*newdm);
1550:   return(0);
1551: }


1554: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1555: {

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

1562:   ((DM_Moab*)dm->data)->bs = 1;
1563:   ((DM_Moab*)dm->data)->numFields = 1;
1564:   ((DM_Moab*)dm->data)->n = 0;
1565:   ((DM_Moab*)dm->data)->nloc = 0;
1566:   ((DM_Moab*)dm->data)->nghost = 0;
1567:   ((DM_Moab*)dm->data)->nele = 0;
1568:   ((DM_Moab*)dm->data)->neleloc = 0;
1569:   ((DM_Moab*)dm->data)->neleghost = 0;
1570:   ((DM_Moab*)dm->data)->ltog_map = NULL;
1571:   ((DM_Moab*)dm->data)->ltog_sendrecv = NULL;

1573:   ((DM_Moab*)dm->data)->refct = 1;
1574:   ((DM_Moab*)dm->data)->parent = NULL;
1575:   ((DM_Moab*)dm->data)->vlocal = new moab::Range();
1576:   ((DM_Moab*)dm->data)->vowned = new moab::Range();
1577:   ((DM_Moab*)dm->data)->vghost = new moab::Range();
1578:   ((DM_Moab*)dm->data)->elocal = new moab::Range();
1579:   ((DM_Moab*)dm->data)->eghost = new moab::Range();

1581:   DMInitialize_Moab(dm);
1582:   return(0);
1583: }