Actual source code: dmmoab.cxx

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

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

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

 17:   Reference: https://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html

 19:   Level: intermediate

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

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

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

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

 55:   Collective

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

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

 63:   Level: beginner

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

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

 79:   Collective

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

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

 91:   Level: intermediate

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

101:   PetscFunctionBegin;
102:   PetscAssertPointer(dmb, 6);

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

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

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

120: #ifdef MOAB_HAVE_MPI
121:   moab::EntityHandle partnset;

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

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

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

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

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

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

162: #ifdef MOAB_HAVE_MPI

164: /*@C
165:   DMMoabGetParallelComm - Get the ParallelComm used with this DMMoab

167:   Collective

169:   Input Parameter:
170: . dm    - The DMMoab object being set

172:   Output Parameter:
173: . pcomm - The ParallelComm for the DMMoab

175:   Level: beginner

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

187: #endif /* MOAB_HAVE_MPI */

189: /*@C
190:   DMMoabSetInterface - Set the MOAB instance used with this DMMoab

192:   Collective

194:   Input Parameters:
195: + dm      - The DMMoab object being set
196: - mbiface - The MOAB instance being set on this DMMoab

198:   Level: beginner

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

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

217: /*@C
218:   DMMoabGetInterface - Get the MOAB instance used with this DMMoab

220:   Collective

222:   Input Parameter:
223: . dm - The DMMoab object being set

225:   Output Parameter:
226: . mbiface - The MOAB instance set on this DMMoab

228:   Level: beginner

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

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

245: /*@C
246:   DMMoabSetLocalVertices - Set the entities having DOFs on this DMMoab

248:   Collective

250:   Input Parameters:
251: + dm    - The DMMoab object being set
252: - range - The entities treated by this DMMoab

254:   Level: beginner

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

263:   PetscFunctionBegin;
265:   dmmoab->vlocal->clear();
266:   dmmoab->vowned->clear();

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

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

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

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

297: /*@C
298:   DMMoabGetAllVertices - Get the entities having DOFs on this DMMoab

300:   Collective

302:   Input Parameter:
303: . dm - The DMMoab object being set

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

308:   Level: beginner

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

320: /*@C
321:   DMMoabGetLocalVertices - Get the entities having DOFs on this DMMoab

323:   Collective

325:   Input Parameter:
326: . dm - The DMMoab object being set

328:   Output Parameters:
329: + owned - The owned vertex entities in this DMMoab
330: - ghost - The ghosted entities (non-owned) stored locally in this partition

332:   Level: beginner

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

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

348:   Collective

350:   Input Parameter:
351: . dm - The DMMoab object being set

353:   Output Parameter:
354: . range - The entities owned locally

356:   Level: beginner

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

368: /*@C
369:   DMMoabSetLocalElements - Set the entities having DOFs on this DMMoab

371:   Collective

373:   Input Parameters:
374: + dm    - The DMMoab object being set
375: - range - The entities treated by this DMMoab

377:   Level: beginner

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

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

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

410:   Collective

412:   Input Parameters:
413: + dm      - The DMMoab object being set
414: - ltogtag - The MOAB tag used for local to global ids

416:   Level: beginner

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

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

431:   Collective

433:   Input Parameter:
434: . dm - The DMMoab object being set

436:   Output Parameter:
437: . ltog_tag - The MOAB tag used for local to global ids

439:   Level: beginner

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

451: /*@C
452:   DMMoabSetBlockSize - Set the block size used with this DMMoab

454:   Collective

456:   Input Parameters:
457: + dm - The DMMoab object being set
458: - bs - The block size used with this DMMoab

460:   Level: beginner

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

472: /*@C
473:   DMMoabGetBlockSize - Get the block size used with this DMMoab

475:   Collective

477:   Input Parameter:
478: . dm - The DMMoab object being set

480:   Output Parameter:
481: . bs - The block size used with this DMMoab

483:   Level: beginner

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

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

498:   Collective on dm

500:   Input Parameter:
501: . dm - The DMMoab object being set

503:   Output Parameters:
504: + neg - The number of global elements in the DMMoab instance
505: - nvg - The number of global vertices in the DMMoab instance

507:   Level: beginner

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

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

523:   Collective on dm

525:   Input Parameter:
526: . dm - The DMMoab object being set

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

534:   Level: beginner

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

549: /*@C
550:   DMMoabGetOffset - Get the local offset for the global vector

552:   Collective

554:   Input Parameter:
555: . dm - The DMMoab object being set

557:   Output Parameter:
558: . offset - The local offset for the global vector

560:   Level: beginner

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

572: /*@C
573:   DMMoabGetDimension - Get the dimension of the DM Mesh

575:   Collective

577:   Input Parameter:
578: . dm - The DMMoab object

580:   Output Parameter:
581: . dim - The dimension of DM

583:   Level: beginner

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

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

599:   Collective on dm

601:   Input Parameter:
602: . dm - The DMMoab object being set

604:   Output Parameter:
605: . nlevel - The current mesh hierarchy level

607:   Level: beginner

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

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

622:   Collective

624:   Input Parameters:
625: + dm      - The DMMoab object
626: - ehandle - The element entity handle

628:   Output Parameter:
629: . mat - The material ID for the current entity

631:   Level: beginner

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

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

648: /*@C
649:   DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities

651:   Collective

653:   Input Parameters:
654: + dm    - The DMMoab object
655: . nconn - Number of entities whose coordinates are needed
656: - conn  - The vertex entity handles

658:   Output Parameter:
659: . vpos - The coordinates of the requested vertex entities

661:   Level: beginner

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

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

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

687: /*@C
688:   DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity

690:   Collective

692:   Input Parameters:
693: + dm      - The DMMoab object
694: - vhandle - Vertex entity handle

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

700:   Level: beginner

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

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

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

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

729: /*@C
730:   DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity

732:   Collective

734:   Input Parameters:
735: + dm      - The DMMoab object
736: . ehandle - Vertex entity handle
737: . nconn   - Number of entities whose coordinates are needed
738: - conn    - The vertex entity handles

740:   Level: beginner

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

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

755: /*@C
756:   DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity

758:   Collective

760:   Input Parameters:
761: + dm      - The DMMoab object
762: - ehandle - Vertex entity handle

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

768:   Level: beginner

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

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

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

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

796:   Collective

798:   Input Parameters:
799: + dm  - The DMMoab object
800: - ent - Entity handle

802:   Output Parameter:
803: . ent_on_boundary - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise

805:   Level: beginner

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

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

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

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

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

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

843:   Input Parameters:
844: + dm    - The DMMoab object
845: . nconn - Number of handles
846: - cnt   - Array of entity handles

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

851:   Level: beginner

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

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

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

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

873:   Input Parameter:
874: . dm - The DMMoab object

876:   Output Parameters:
877: + bdvtx   - Boundary vertices
878: . bdelems - Boundary elements
879: - bdfaces - Boundary faces

881:   Level: beginner

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

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

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

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

905:   PetscFunctionBegin;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1049:     dmmoab->neleloc   = dmmoab->elocal->size();
1050:     dmmoab->neleghost = dmmoab->eghost->size();

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1248:   PetscFunctionReturn(PETSC_SUCCESS);
1249: }

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

1254:   Collective

1256:   Input Parameters:
1257: + dm     - The `DM` object
1258: . coords - The connectivity of the element
1259: - nverts - The number of vertices that form the element

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

1264:   Level: beginner

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

1274:   PetscFunctionBegin;
1276:   PetscAssertPointer(coords, 2);

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

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

1286:   if (overts) *overts = verts;
1287:   PetscFunctionReturn(PETSC_SUCCESS);
1288: }

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

1293:   Collective

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

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

1304:   Level: beginner

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

1314:   PetscFunctionBegin;
1316:   PetscAssertPointer(conn, 3);

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

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

1326:   if (oelem) *oelem = elem;
1327:   PetscFunctionReturn(PETSC_SUCCESS);
1328: }

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

1335:   Collective

1337:   Input Parameters:
1338: . dm - The DM object

1340:   Output Parameter:
1341: . newdm - The sub DM object with updated set information

1343:   Level: advanced

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

1353:   PetscFunctionBegin;

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

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

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

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

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

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

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

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

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

1421: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1422: {
1423:   PetscFunctionReturn(PETSC_SUCCESS);
1424: }

1426: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1427: {
1428:   PetscFunctionReturn(PETSC_SUCCESS);
1429: }

1431: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1432: {
1433:   PetscBool iascii, ishdf5, isvtk;

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

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

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

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

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

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

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

1528:   PetscCall(DMInitialize_Moab(dm));
1529:   PetscFunctionReturn(PETSC_SUCCESS);
1530: }