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 {cite}`moabwebsite`.
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: Level: intermediate
19: .seealso: `DMMOAB`, `DMType`, `DMMoabCreate()`, `DMCreate()`, `DMSetType()`, `DMMoabCreateMoab()`
20: M*/
22: /* External function declarations here */
23: PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
24: PETSC_EXTERN PetscErrorCode DMCreateDefaultConstraints_Moab(DM dm);
25: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
26: PETSC_EXTERN PetscErrorCode DMCreateCoordinateDM_Moab(DM dm, DM *cdm);
27: PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM *dmRefined);
28: PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM *dmCoarsened);
29: PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmRefined[]);
30: PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmCoarsened[]);
31: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm);
32: PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM, Vec *);
33: PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM, Vec *);
34: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
35: PETSC_EXTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM, Vec, InsertMode, Vec);
36: PETSC_EXTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM, Vec, InsertMode, Vec);
37: PETSC_EXTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM, Vec, InsertMode, Vec);
38: PETSC_EXTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM, Vec, InsertMode, Vec);
40: /* Un-implemented routines */
41: /*
42: PETSC_EXTERN PetscErrorCode DMCreatelocalsection_Moab(DM dm);
43: PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dmCoarse, DM dmFine, Mat *mat);
44: PETSC_EXTERN PetscErrorCode DMLoad_Moab(DM dm, PetscViewer viewer);
45: PETSC_EXTERN PetscErrorCode DMGetDimPoints_Moab(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd);
46: PETSC_EXTERN PetscErrorCode DMCreateSubDM_Moab(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
47: PETSC_EXTERN PetscErrorCode DMLocatePoints_Moab(DM dm, Vec v, IS *cellIS);
48: */
50: /*@C
51: DMMoabCreate - Creates a `DMMOAB` object, which encapsulates a moab instance
53: Collective
55: Input Parameter:
56: . comm - The communicator for the `DMMOAB` object
58: Output Parameter:
59: . dmb - The `DMMOAB` object
61: Level: beginner
63: .seealso: `DMMOAB`, `DMMoabCreateMoab()`
64: @*/
65: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
66: {
67: PetscFunctionBegin;
68: PetscAssertPointer(dmb, 2);
69: PetscCall(DMCreate(comm, dmb));
70: PetscCall(DMSetType(*dmb, DMMOAB));
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: /*@C
75: DMMoabCreateMoab - Creates a `DMMOAB` object, optionally from an instance and other data
77: Collective
79: Input Parameters:
80: + comm - The communicator for the `DMMOAB` object
81: . mbiface - (ptr to) the `DMMOAB` Instance; if passed in `NULL`, MOAB instance is created inside PETSc, and destroyed
82: along with the `DMMOAB`
83: . ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
84: - range - If non-`NULL`, contains range of entities to which DOFs will be assigned
86: Output Parameter:
87: . dmb - The `DMMOAB` object
89: Level: intermediate
91: .seealso: `DMMOAB`, `DMMoabCreate()`
92: @*/
93: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
94: {
95: moab::ErrorCode merr;
96: DM dmmb;
97: DM_Moab *dmmoab;
99: PetscFunctionBegin;
100: PetscAssertPointer(dmb, 6);
102: PetscCall(DMMoabCreate(comm, &dmmb));
103: dmmoab = (DM_Moab *)dmmb->data;
105: if (!mbiface) {
106: dmmoab->mbiface = new moab::Core();
107: dmmoab->icreatedinstance = PETSC_TRUE;
108: } else {
109: dmmoab->mbiface = mbiface;
110: dmmoab->icreatedinstance = PETSC_FALSE;
111: }
113: /* by default the fileset = root set. This set stores the hierarchy of entities belonging to current DM */
114: dmmoab->fileset = 0;
115: dmmoab->hlevel = 0;
116: dmmoab->nghostrings = 0;
118: #ifdef MOAB_HAVE_MPI
119: moab::EntityHandle partnset;
121: /* Create root sets for each mesh. Then pass these
122: to the load_file functions to be populated. */
123: merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, partnset);
124: MBERR("Creating partition set failed", merr);
126: /* Create the parallel communicator object with the partition handle associated with MOAB */
127: dmmoab->pcomm = moab::ParallelComm::get_pcomm(dmmoab->mbiface, partnset, &comm);
128: #endif
130: /* do the remaining initializations for DMMoab */
131: dmmoab->bs = 1;
132: dmmoab->numFields = 1;
133: PetscCall(PetscMalloc(dmmoab->numFields * sizeof(char *), &dmmoab->fieldNames));
134: PetscCall(PetscStrallocpy("DEFAULT", (char **)&dmmoab->fieldNames[0]));
135: dmmoab->rw_dbglevel = 0;
136: dmmoab->partition_by_rank = PETSC_FALSE;
137: dmmoab->extra_read_options[0] = '\0';
138: dmmoab->extra_write_options[0] = '\0';
139: dmmoab->read_mode = READ_PART;
140: dmmoab->write_mode = WRITE_PART;
142: /* set global ID tag handle */
143: if (ltog_tag && *ltog_tag) PetscCall(DMMoabSetLocalToGlobalTag(dmmb, *ltog_tag));
144: else {
145: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
146: MBERRNM(merr);
147: if (ltog_tag) *ltog_tag = dmmoab->ltog_tag;
148: }
150: merr = dmmoab->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dmmoab->material_tag);
151: MBERRNM(merr);
153: /* set the local range of entities (vertices) of interest */
154: if (range) PetscCall(DMMoabSetLocalVertices(dmmb, range));
155: *dmb = dmmb;
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: #ifdef MOAB_HAVE_MPI
161: /*@C
162: DMMoabGetParallelComm - Get the ParallelComm used with this `DMMOAB`
164: Collective
166: Input Parameter:
167: . dm - The `DMMOAB` object being set
169: Output Parameter:
170: . pcomm - The ParallelComm for the `DMMOAB`
172: Level: beginner
174: .seealso: `DMMOAB`, `DMMoabSetInterface()`
175: @*/
176: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
177: {
178: PetscFunctionBegin;
180: *pcomm = ((DM_Moab *)dm->data)->pcomm;
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: #endif /* MOAB_HAVE_MPI */
186: /*@C
187: DMMoabSetInterface - Set the MOAB instance used with this `DMMOAB`
189: Collective
191: Input Parameters:
192: + dm - The `DMMOAB` object being set
193: - mbiface - The MOAB instance being set on this `DMMOAB`
195: Level: beginner
197: .seealso: `DMMOAB`, `DMMoabGetInterface()`
198: @*/
199: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
200: {
201: DM_Moab *dmmoab = (DM_Moab *)dm->data;
203: PetscFunctionBegin;
205: PetscAssertPointer(mbiface, 2);
206: #ifdef MOAB_HAVE_MPI
207: dmmoab->pcomm = NULL;
208: #endif
209: dmmoab->mbiface = mbiface;
210: dmmoab->icreatedinstance = PETSC_FALSE;
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: /*@C
215: DMMoabGetInterface - Get the MOAB instance used with this `DMMOAB`
217: Collective
219: Input Parameter:
220: . dm - The `DMMOAB` object being set
222: Output Parameter:
223: . mbiface - The MOAB instance set on this `DMMOAB`
225: Level: beginner
227: .seealso: `DMMOAB`, `DMMoabSetInterface()`
228: @*/
229: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
230: {
231: static PetscBool cite = PETSC_FALSE;
233: PetscFunctionBegin;
235: 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, "
236: "K. and Stimpson, C. and Ernst, C.},\n year = {2004}, note = {Report}\n}\n",
237: &cite));
238: *mbiface = ((DM_Moab *)dm->data)->mbiface;
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: /*@C
243: DMMoabSetLocalVertices - Set the entities having DOFs on this `DMMOAB`
245: Collective
247: Input Parameters:
248: + dm - The `DMMOAB` object being set
249: - range - The entities treated by this `DMMOAB`
251: Level: beginner
253: .seealso: `DMMOAB`, `DMMoabGetAllVertices()`
254: @*/
255: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
256: {
257: moab::Range tmpvtxs;
258: DM_Moab *dmmoab = (DM_Moab *)dm->data;
260: PetscFunctionBegin;
262: dmmoab->vlocal->clear();
263: dmmoab->vowned->clear();
265: dmmoab->vlocal->insert(range->begin(), range->end());
267: #ifdef MOAB_HAVE_MPI
268: moab::ErrorCode merr;
269: /* filter based on parallel status */
270: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
271: MBERRNM(merr);
273: /* filter all the non-owned and shared entities out of the list */
274: tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
275: merr = dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
276: MBERRNM(merr);
277: tmpvtxs = moab::subtract(tmpvtxs, *dmmoab->vghost);
278: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
279: #else
280: *dmmoab->vowned = *dmmoab->vlocal;
281: #endif
283: /* compute and cache the sizes of local and ghosted entities */
284: dmmoab->nloc = dmmoab->vowned->size();
285: dmmoab->nghost = dmmoab->vghost->size();
286: #ifdef MOAB_HAVE_MPI
287: PetscCallMPI(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
288: #else
289: dmmoab->n = dmmoab->nloc;
290: #endif
291: PetscFunctionReturn(PETSC_SUCCESS);
292: }
294: /*@C
295: DMMoabGetAllVertices - Get the entities having DOFs on this `DMMOAB`
297: Collective
299: Input Parameter:
300: . dm - The `DMMOAB` object being set
302: Output Parameter:
303: . local - The local vertex entities in this `DMMOAB` = (owned+ghosted)
305: Level: beginner
307: .seealso: `DMMOAB`, `DMMoabGetLocalVertices()`
308: @*/
309: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
310: {
311: PetscFunctionBegin;
313: if (local) *local = *((DM_Moab *)dm->data)->vlocal;
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: /*@C
318: DMMoabGetLocalVertices - Get the entities having DOFs on this `DMMOAB`
320: Collective
322: Input Parameter:
323: . dm - The `DMMOAB` object being set
325: Output Parameters:
326: + owned - The owned vertex entities in this `DMMOAB`
327: - ghost - The ghosted entities (non-owned) stored locally in this partition
329: Level: beginner
331: .seealso: `DMMOAB`, `DMMoabGetAllVertices()`
332: @*/
333: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
334: {
335: PetscFunctionBegin;
337: if (owned) *owned = ((DM_Moab *)dm->data)->vowned;
338: if (ghost) *ghost = ((DM_Moab *)dm->data)->vghost;
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: /*@C
343: DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned
345: Collective
347: Input Parameter:
348: . dm - The `DMMOAB` object being set
350: Output Parameter:
351: . range - The entities owned locally
353: Level: beginner
355: .seealso: `DMMOAB`, `DMMoabSetLocalElements()`
356: @*/
357: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
358: {
359: PetscFunctionBegin;
361: if (range) *range = ((DM_Moab *)dm->data)->elocal;
362: PetscFunctionReturn(PETSC_SUCCESS);
363: }
365: /*@C
366: DMMoabSetLocalElements - Set the entities having DOFs on this `DMMOAB`
368: Collective
370: Input Parameters:
371: + dm - The `DMMOAB` object being set
372: - range - The entities treated by this `DMMOAB`
374: Level: beginner
376: .seealso: `DMMOAB`, `DMMoabGetLocalElements()`
377: @*/
378: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
379: {
380: DM_Moab *dmmoab = (DM_Moab *)dm->data;
382: PetscFunctionBegin;
384: dmmoab->elocal->clear();
385: dmmoab->eghost->clear();
386: dmmoab->elocal->insert(range->begin(), range->end());
387: #ifdef MOAB_HAVE_MPI
388: moab::ErrorCode merr;
389: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
390: MBERRNM(merr);
391: *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
392: #endif
393: dmmoab->neleloc = dmmoab->elocal->size();
394: dmmoab->neleghost = dmmoab->eghost->size();
395: #ifdef MOAB_HAVE_MPI
396: PetscCallMPI(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
397: PetscCall(PetscInfo(dm, "Created %" PetscInt_FMT " local and %" PetscInt_FMT " global elements.\n", dmmoab->neleloc, dmmoab->nele));
398: #else
399: dmmoab->nele = dmmoab->neleloc;
400: #endif
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
404: /*@C
405: DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
407: Collective
409: Input Parameters:
410: + dm - The `DMMOAB` object being set
411: - ltogtag - The `DMMOAB` tag used for local to global ids
413: Level: beginner
415: .seealso: `DMMOAB`, `DMMoabGetLocalToGlobalTag()`
416: @*/
417: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
418: {
419: PetscFunctionBegin;
421: ((DM_Moab *)dm->data)->ltog_tag = ltogtag;
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: /*@C
426: DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
428: Collective
430: Input Parameter:
431: . dm - The `DMMOAB` object being set
433: Output Parameter:
434: . ltog_tag - The MOAB tag used for local to global ids
436: Level: beginner
438: .seealso: `DMMOAB`, `DMMoabSetLocalToGlobalTag()`
439: @*/
440: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
441: {
442: PetscFunctionBegin;
444: *ltog_tag = ((DM_Moab *)dm->data)->ltog_tag;
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: /*@C
449: DMMoabSetBlockSize - Set the block size used with this `DMMOAB`
451: Collective
453: Input Parameters:
454: + dm - The `DMMOAB` object being set
455: - bs - The block size used with this `DMMOAB`
457: Level: beginner
459: .seealso: `DMMOAB`, `DMMoabGetBlockSize()`
460: @*/
461: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
462: {
463: PetscFunctionBegin;
465: ((DM_Moab *)dm->data)->bs = bs;
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: /*@C
470: DMMoabGetBlockSize - Get the block size used with this `DMMOAB`
472: Collective
474: Input Parameter:
475: . dm - The `DMMOAB` object being set
477: Output Parameter:
478: . bs - The block size used with this `DMMOAB`
480: Level: beginner
482: .seealso: `DMMOAB`, `DMMoabSetBlockSize()`
483: @*/
484: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
485: {
486: PetscFunctionBegin;
488: *bs = ((DM_Moab *)dm->data)->bs;
489: PetscFunctionReturn(PETSC_SUCCESS);
490: }
492: /*@C
493: DMMoabGetSize - Get the global vertex size used with this `DMMOAB`
495: Collective
497: Input Parameter:
498: . dm - The `DMMOAB` object being set
500: Output Parameters:
501: + neg - The number of global elements in the `DMMOAB` instance
502: - nvg - The number of global vertices in the `DMMOAB` instance
504: Level: beginner
506: .seealso: `DMMOAB`, `DMMoabGetLocalSize()`
507: @*/
508: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
509: {
510: PetscFunctionBegin;
512: if (neg) *neg = ((DM_Moab *)dm->data)->nele;
513: if (nvg) *nvg = ((DM_Moab *)dm->data)->n;
514: PetscFunctionReturn(PETSC_SUCCESS);
515: }
517: /*@C
518: DMMoabGetLocalSize - Get the local and ghosted vertex size used with this `DMMOAB`
520: Collective
522: Input Parameter:
523: . dm - The `DMMOAB` object being set
525: Output Parameters:
526: + nel - The number of owned elements in this processor
527: . neg - The number of ghosted elements in this processor
528: . nvl - The number of owned vertices in this processor
529: - nvg - The number of ghosted vertices in this processor
531: Level: beginner
533: .seealso: `DMMOAB`, `DMMoabGetSize()`
534: @*/
535: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
536: {
537: PetscFunctionBegin;
539: if (nel) *nel = ((DM_Moab *)dm->data)->neleloc;
540: if (neg) *neg = ((DM_Moab *)dm->data)->neleghost;
541: if (nvl) *nvl = ((DM_Moab *)dm->data)->nloc;
542: if (nvg) *nvg = ((DM_Moab *)dm->data)->nghost;
543: PetscFunctionReturn(PETSC_SUCCESS);
544: }
546: /*@C
547: DMMoabGetOffset - Get the local offset for the global vector
549: Collective
551: Input Parameter:
552: . dm - The `DMMOAB` object being set
554: Output Parameter:
555: . offset - The local offset for the global vector
557: Level: beginner
559: .seealso: `DMMOAB`, `DMMoabGetDimension()`
560: @*/
561: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
562: {
563: PetscFunctionBegin;
565: *offset = ((DM_Moab *)dm->data)->vstart;
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: /*@C
570: DMMoabGetDimension - Get the dimension of the `DM` Mesh
572: Collective
574: Input Parameter:
575: . dm - The `DMMOAB` object
577: Output Parameter:
578: . dim - The dimension of `DM`
580: Level: beginner
582: .seealso: `DMMOAB`, `DMMoabGetOffset()`
583: @*/
584: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
585: {
586: PetscFunctionBegin;
588: *dim = ((DM_Moab *)dm->data)->dim;
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: /*@C
593: DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
594: generated through uniform refinement.
596: Collective
598: Input Parameter:
599: . dm - The `DMMOAB` object being set
601: Output Parameter:
602: . nlevel - The current mesh hierarchy level
604: Level: beginner
606: .seealso: `DMMOAB`, `DMMoabGetMaterialBlock()`
607: @*/
608: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
609: {
610: PetscFunctionBegin;
612: if (nlevel) *nlevel = ((DM_Moab *)dm->data)->hlevel;
613: PetscFunctionReturn(PETSC_SUCCESS);
614: }
616: /*@C
617: DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh
619: Collective
621: Input Parameters:
622: + dm - The `DMMOAB` object
623: - ehandle - The element entity handle
625: Output Parameter:
626: . mat - The material ID for the current entity
628: Level: beginner
630: .seealso: `DMMOAB`, `DMMoabGetHierarchyLevel()`
631: @*/
632: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
633: {
634: DM_Moab *dmmoab;
636: PetscFunctionBegin;
638: if (*mat) {
639: dmmoab = (DM_Moab *)dm->data;
640: *mat = dmmoab->materials[dmmoab->elocal->index(ehandle)];
641: }
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: /*@C
646: DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities
648: Collective
650: Input Parameters:
651: + dm - The `DMMOAB` object
652: . nconn - Number of entities whose coordinates are needed
653: - conn - The vertex entity handles
655: Output Parameter:
656: . vpos - The coordinates of the requested vertex entities
658: Level: beginner
660: .seealso: `DMMOAB`, `DMMoabGetVertexConnectivity()`
661: @*/
662: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
663: {
664: DM_Moab *dmmoab;
665: moab::ErrorCode merr;
667: PetscFunctionBegin;
669: PetscAssertPointer(conn, 3);
670: PetscAssertPointer(vpos, 4);
671: dmmoab = (DM_Moab *)dm->data;
673: /* Get connectivity information in MOAB canonical ordering */
674: if (dmmoab->hlevel) {
675: merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle *>(conn), nconn, dmmoab->hlevel, vpos);
676: MBERRNM(merr);
677: } else {
678: merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);
679: MBERRNM(merr);
680: }
681: PetscFunctionReturn(PETSC_SUCCESS);
682: }
684: /*@C
685: DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity
687: Collective
689: Input Parameters:
690: + dm - The `DMMOAB` object
691: - vhandle - Vertex entity handle
693: Output Parameters:
694: + nconn - Number of entities whose coordinates are needed
695: - conn - The vertex entity handles
697: Level: beginner
699: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabRestoreVertexConnectivity()`
700: @*/
701: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt *nconn, moab::EntityHandle **conn)
702: {
703: DM_Moab *dmmoab;
704: std::vector<moab::EntityHandle> adj_entities, connect;
705: moab::ErrorCode merr;
707: PetscFunctionBegin;
709: PetscAssertPointer(conn, 4);
710: dmmoab = (DM_Moab *)dm->data;
712: /* Get connectivity information in MOAB canonical ordering */
713: merr = dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION);
714: MBERRNM(merr);
715: merr = dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect);
716: MBERRNM(merr);
718: if (conn) {
719: PetscCall(PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn));
720: PetscCall(PetscArraycpy(*conn, &connect[0], connect.size()));
721: }
722: if (nconn) *nconn = connect.size();
723: PetscFunctionReturn(PETSC_SUCCESS);
724: }
726: /*@C
727: DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity
729: Collective
731: Input Parameters:
732: + dm - The `DMMOAB` object
733: . ehandle - Vertex entity handle
734: . nconn - Number of entities whose coordinates are needed
735: - conn - The vertex entity handles
737: Level: beginner
739: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`
740: @*/
741: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, moab::EntityHandle **conn)
742: {
743: PetscFunctionBegin;
745: PetscAssertPointer(conn, 4);
747: if (conn) PetscCall(PetscFree(*conn));
748: if (nconn) *nconn = 0;
749: PetscFunctionReturn(PETSC_SUCCESS);
750: }
752: /*@C
753: DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity
755: Collective
757: Input Parameters:
758: + dm - The `DMMOAB` object
759: - ehandle - Vertex entity handle
761: Output Parameters:
762: + nconn - Number of entities whose coordinates are needed
763: - conn - The vertex entity handles
765: Level: beginner
767: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`, `DMMoabRestoreVertexConnectivity()`
768: @*/
769: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, const moab::EntityHandle **conn)
770: {
771: DM_Moab *dmmoab;
772: const moab::EntityHandle *connect;
773: std::vector<moab::EntityHandle> vconn;
774: moab::ErrorCode merr;
775: PetscInt nnodes;
777: PetscFunctionBegin;
779: PetscAssertPointer(conn, 4);
780: dmmoab = (DM_Moab *)dm->data;
782: /* Get connectivity information in MOAB canonical ordering */
783: merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes);
784: MBERRNM(merr);
785: if (conn) *conn = connect;
786: if (nconn) *nconn = nnodes;
787: PetscFunctionReturn(PETSC_SUCCESS);
788: }
790: /*@C
791: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
793: Collective
795: Input Parameters:
796: + dm - The `DMMOAB` object
797: - ent - Entity handle
799: Output Parameter:
800: . ent_on_boundary - `PETSC_TRUE` if entity on boundary; `PETSC_FALSE` otherwise
802: Level: beginner
804: .seealso: `DMMOAB`, `DMMoabCheckBoundaryVertices()`
805: @*/
806: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool *ent_on_boundary)
807: {
808: moab::EntityType etype;
809: DM_Moab *dmmoab;
810: PetscInt edim;
812: PetscFunctionBegin;
814: PetscAssertPointer(ent_on_boundary, 3);
815: dmmoab = (DM_Moab *)dm->data;
817: /* get the entity type and handle accordingly */
818: etype = dmmoab->mbiface->type_from_handle(ent);
819: PetscCheck(etype < moab::MBPOLYHEDRON, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Entity type on the boundary skin is invalid. EntityType = %" PetscInt_FMT, etype);
821: /* get the entity dimension */
822: edim = dmmoab->mbiface->dimension_from_handle(ent);
824: *ent_on_boundary = PETSC_FALSE;
825: if (etype == moab::MBVERTEX && edim == 0) {
826: *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
827: } else {
828: if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
829: if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
830: } else { /* next check the lower-dimensional faces */
831: if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
832: }
833: }
834: PetscFunctionReturn(PETSC_SUCCESS);
835: }
837: /*@C
838: DMMoabCheckBoundaryVertices - Check whether a given entity is on the boundary (vertex, edge, face, element)
840: Input Parameters:
841: + dm - The `DMMOAB` object
842: . nconn - Number of handles
843: - cnt - Array of entity handles
845: Output Parameter:
846: . isbdvtx - Array of boundary markers - `PETSC_TRUE` if entity on boundary; `PETSC_FALSE` otherwise
848: Level: beginner
850: .seealso: `DMMOAB`, `DMMoabIsEntityOnBoundary()`
851: @*/
852: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool *isbdvtx)
853: {
854: DM_Moab *dmmoab;
855: PetscInt i;
857: PetscFunctionBegin;
859: PetscAssertPointer(cnt, 3);
860: PetscAssertPointer(isbdvtx, 4);
861: dmmoab = (DM_Moab *)dm->data;
863: for (i = 0; i < nconn; ++i) isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
864: PetscFunctionReturn(PETSC_SUCCESS);
865: }
867: /*@C
868: DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary
870: Input Parameter:
871: . dm - The `DMMOAB` object
873: Output Parameters:
874: + bdvtx - Boundary vertices
875: . bdelems - Boundary elements
876: - bdfaces - Boundary faces
878: Level: beginner
880: .seealso: `DMMOAB`, `DMMoabCheckBoundaryVertices()`, `DMMoabIsEntityOnBoundary()`
881: @*/
882: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range **bdelems, const moab::Range **bdfaces)
883: {
884: DM_Moab *dmmoab;
886: PetscFunctionBegin;
888: dmmoab = (DM_Moab *)dm->data;
890: if (bdvtx) *bdvtx = dmmoab->bndyvtx;
891: if (bdfaces) *bdfaces = dmmoab->bndyfaces;
892: if (bdelems) *bdfaces = dmmoab->bndyelems;
893: PetscFunctionReturn(PETSC_SUCCESS);
894: }
896: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
897: {
898: PetscInt i;
899: moab::ErrorCode merr;
900: DM_Moab *dmmoab = (DM_Moab *)dm->data;
902: PetscFunctionBegin;
905: dmmoab->refct--;
906: if (!dmmoab->refct) {
907: delete dmmoab->vlocal;
908: delete dmmoab->vowned;
909: delete dmmoab->vghost;
910: delete dmmoab->elocal;
911: delete dmmoab->eghost;
912: delete dmmoab->bndyvtx;
913: delete dmmoab->bndyfaces;
914: delete dmmoab->bndyelems;
916: PetscCall(PetscFree(dmmoab->gsindices));
917: PetscCall(PetscFree2(dmmoab->gidmap, dmmoab->lidmap));
918: PetscCall(PetscFree(dmmoab->dfill));
919: PetscCall(PetscFree(dmmoab->ofill));
920: PetscCall(PetscFree(dmmoab->materials));
921: if (dmmoab->fieldNames) {
922: for (i = 0; i < dmmoab->numFields; i++) PetscCall(PetscFree(dmmoab->fieldNames[i]));
923: PetscCall(PetscFree(dmmoab->fieldNames));
924: }
926: if (dmmoab->nhlevels) {
927: PetscCall(PetscFree(dmmoab->hsets));
928: dmmoab->nhlevels = 0;
929: if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
930: dmmoab->hierarchy = NULL;
931: }
933: if (dmmoab->icreatedinstance) {
934: delete dmmoab->pcomm;
935: merr = dmmoab->mbiface->delete_mesh();
936: MBERRNM(merr);
937: delete dmmoab->mbiface;
938: }
939: dmmoab->mbiface = NULL;
940: #ifdef MOAB_HAVE_MPI
941: dmmoab->pcomm = NULL;
942: #endif
943: PetscCall(VecScatterDestroy(&dmmoab->ltog_sendrecv));
944: PetscCall(ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map));
945: PetscCall(PetscFree(dm->data));
946: }
947: PetscFunctionReturn(PETSC_SUCCESS);
948: }
950: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(DM dm, PetscOptionItems PetscOptionsObject)
951: {
952: DM_Moab *dmmoab = (DM_Moab *)dm->data;
954: PetscFunctionBegin;
955: PetscOptionsHeadBegin(PetscOptionsObject, "DMMoab Options");
956: PetscCall(PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL, 0));
957: 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));
958: /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
959: 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));
960: 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));
961: PetscCall(PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum *)&dmmoab->read_mode, NULL));
962: PetscCall(PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum *)&dmmoab->write_mode, NULL));
963: PetscOptionsHeadEnd();
964: PetscFunctionReturn(PETSC_SUCCESS);
965: }
967: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
968: {
969: moab::ErrorCode merr;
970: Vec local, global;
971: IS from, to;
972: moab::Range::iterator iter;
973: PetscInt i, j, f, bs, vent, totsize, *lgmap;
974: DM_Moab *dmmoab = (DM_Moab *)dm->data;
975: moab::Range adjs;
977: PetscFunctionBegin;
979: /* Get the local and shared vertices and cache it */
980: PetscCheck(dmmoab->mbiface != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB Interface before calling SetUp.");
981: #ifdef MOAB_HAVE_MPI
982: PetscCheck(dmmoab->pcomm != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB ParallelComm object before calling SetUp.");
983: #endif
985: /* Get the entities recursively in the current part of the mesh, if user did not set the local vertices explicitly */
986: if (dmmoab->vlocal->empty()) {
987: //merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,*dmmoab->vlocal,true);MBERRNM(merr);
988: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, 0, *dmmoab->vlocal, false);
989: MBERRNM(merr);
991: #ifdef MOAB_HAVE_MPI
992: /* filter based on parallel status */
993: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
994: MBERRNM(merr);
996: /* filter all the non-owned and shared entities out of the list */
997: // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
998: adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
999: merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
1000: MBERRNM(merr);
1001: adjs = moab::subtract(adjs, *dmmoab->vghost);
1002: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
1003: #else
1004: *dmmoab->vowned = *dmmoab->vlocal;
1005: #endif
1007: /* compute and cache the sizes of local and ghosted entities */
1008: dmmoab->nloc = dmmoab->vowned->size();
1009: dmmoab->nghost = dmmoab->vghost->size();
1011: #ifdef MOAB_HAVE_MPI
1012: PetscCallMPI(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
1013: 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));
1014: #else
1015: dmmoab->n = dmmoab->nloc;
1016: #endif
1017: }
1019: {
1020: /* get the information about the local elements in the mesh */
1021: dmmoab->eghost->clear();
1023: /* first decipher the leading dimension */
1024: for (i = 3; i > 0; i--) {
1025: dmmoab->elocal->clear();
1026: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false);
1027: MBERRNM(merr);
1029: /* store the current mesh dimension */
1030: if (dmmoab->elocal->size()) {
1031: dmmoab->dim = i;
1032: break;
1033: }
1034: }
1036: PetscCall(DMSetDimension(dm, dmmoab->dim));
1038: #ifdef MOAB_HAVE_MPI
1039: /* filter the ghosted and owned element list */
1040: *dmmoab->eghost = *dmmoab->elocal;
1041: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1042: MBERRNM(merr);
1043: *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
1044: #endif
1046: dmmoab->neleloc = dmmoab->elocal->size();
1047: dmmoab->neleghost = dmmoab->eghost->size();
1049: #ifdef MOAB_HAVE_MPI
1050: PetscCallMPI(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
1051: PetscCall(PetscInfo(NULL, "%d-dim elements: owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost));
1052: #else
1053: dmmoab->nele = dmmoab->neleloc;
1054: #endif
1055: }
1057: bs = dmmoab->bs;
1058: if (!dmmoab->ltog_tag) {
1059: /* Get the global ID tag. The global ID tag is applied to each
1060: vertex. It acts as an global identifier which MOAB uses to
1061: assemble the individual pieces of the mesh */
1062: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
1063: MBERRNM(merr);
1064: }
1066: totsize = dmmoab->vlocal->size();
1067: 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);
1068: PetscCall(PetscCalloc1(totsize, &dmmoab->gsindices));
1069: {
1070: /* first get the local indices */
1071: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]);
1072: MBERRNM(merr);
1073: if (dmmoab->nghost) { /* next get the ghosted indices */
1074: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]);
1075: MBERRNM(merr);
1076: }
1078: /* find out the local and global minima of GLOBAL_ID */
1079: dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1080: for (i = 0; i < totsize; ++i) {
1081: if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1082: if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1083: }
1085: PetscCallMPI(MPIU_Allreduce(&dmmoab->lminmax[0], &dmmoab->gminmax[0], 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm));
1086: PetscCallMPI(MPIU_Allreduce(&dmmoab->lminmax[1], &dmmoab->gminmax[1], 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm));
1088: /* set the GID map */
1089: for (i = 0; i < totsize; ++i) dmmoab->gsindices[i] -= dmmoab->gminmax[0]; /* zero based index needed for IS */
1090: dmmoab->lminmax[0] -= dmmoab->gminmax[0];
1091: dmmoab->lminmax[1] -= dmmoab->gminmax[0];
1093: 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]));
1094: }
1095: 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,
1096: dmmoab->numFields);
1098: {
1099: dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1100: dmmoab->seqend = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1101: PetscCall(PetscInfo(NULL, "SEQUENCE: Local [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "]\n", dmmoab->seqstart, dmmoab->seqend));
1103: PetscCall(PetscMalloc2(dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->gidmap, dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->lidmap));
1104: PetscCall(PetscMalloc1(totsize * dmmoab->numFields, &lgmap));
1106: i = j = 0;
1107: /* set the owned vertex data first */
1108: for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1109: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1110: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1111: dmmoab->lidmap[vent] = i;
1112: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1113: }
1114: /* next arrange all the ghosted data information */
1115: for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1116: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1117: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1118: dmmoab->lidmap[vent] = i;
1119: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1120: }
1122: /* We need to create the Global to Local Vector Scatter Contexts
1123: 1) First create a local and global vector
1124: 2) Create a local and global IS
1125: 3) Create VecScatter and LtoGMapping objects
1126: 4) Cleanup the IS and Vec objects
1127: */
1128: PetscCall(DMCreateGlobalVector(dm, &global));
1129: PetscCall(DMCreateLocalVector(dm, &local));
1131: PetscCall(VecGetOwnershipRange(global, &dmmoab->vstart, &dmmoab->vend));
1133: /* global to local must retrieve ghost points */
1134: PetscCall(ISCreateStride(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, dmmoab->vstart, 1, &from));
1135: PetscCall(ISSetBlockSize(from, bs));
1137: PetscCall(ISCreateGeneral(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, &lgmap[0], PETSC_COPY_VALUES, &to));
1138: PetscCall(ISSetBlockSize(to, bs));
1140: if (!dmmoab->ltog_map) {
1141: /* create to the local to global mapping for vectors in order to use VecSetValuesLocal */
1142: PetscCall(ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm, dmmoab->bs, totsize * dmmoab->numFields, lgmap, PETSC_COPY_VALUES, &dmmoab->ltog_map));
1143: }
1145: /* now create the scatter object from local to global vector */
1146: PetscCall(VecScatterCreate(local, from, global, to, &dmmoab->ltog_sendrecv));
1148: /* clean up IS, Vec */
1149: PetscCall(PetscFree(lgmap));
1150: PetscCall(ISDestroy(&from));
1151: PetscCall(ISDestroy(&to));
1152: PetscCall(VecDestroy(&local));
1153: PetscCall(VecDestroy(&global));
1154: }
1156: dmmoab->bndyvtx = new moab::Range();
1157: dmmoab->bndyfaces = new moab::Range();
1158: dmmoab->bndyelems = new moab::Range();
1159: /* skin the boundary and store nodes */
1160: if (!dmmoab->hlevel) {
1161: /* get the skin vertices of boundary faces for the current partition and then filter
1162: the local, boundary faces, vertices and elements alone via PSTATUS flags;
1163: this should not give us any ghosted boundary, but if user needs such a functionality
1164: it would be easy to add it based on the find_skin query below */
1165: moab::Skinner skinner(dmmoab->mbiface);
1167: /* get the entities on the skin - only the faces */
1168: merr = skinner.find_skin(dmmoab->fileset, *dmmoab->elocal, false, *dmmoab->bndyfaces, NULL, true, true, false);
1169: MBERRNM(merr); // 'false' param indicates we want faces back, not vertices
1171: #ifdef MOAB_HAVE_MPI
1172: /* filter all the non-owned and shared entities out of the list */
1173: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1174: MBERRNM(merr);
1175: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT);
1176: MBERRNM(merr);
1177: #endif
1179: /* get all the nodes via connectivity and the parent elements via adjacency information */
1180: merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false);
1181: MBERRNM(merr);
1182: merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION);
1183: MBERRNM(merr);
1184: } else {
1185: /* Let us query the hierarchy manager and get the results directly for this level */
1186: for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1187: moab::EntityHandle elemHandle = *iter;
1188: if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1189: dmmoab->bndyelems->insert(elemHandle);
1190: /* For this boundary element, query the vertices and add them to the list */
1191: std::vector<moab::EntityHandle> connect;
1192: merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect);
1193: MBERRNM(merr);
1194: for (unsigned iv = 0; iv < connect.size(); ++iv)
1195: if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv])) dmmoab->bndyvtx->insert(connect[iv]);
1196: /* Next, let us query the boundary faces and add them also to the list */
1197: std::vector<moab::EntityHandle> faces;
1198: merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim - 1, faces);
1199: MBERRNM(merr);
1200: for (unsigned ifa = 0; ifa < faces.size(); ++ifa)
1201: if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa])) dmmoab->bndyfaces->insert(faces[ifa]);
1202: }
1203: }
1204: #ifdef MOAB_HAVE_MPI
1205: /* filter all the non-owned and shared entities out of the list */
1206: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1207: MBERRNM(merr);
1208: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1209: MBERRNM(merr);
1210: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1211: MBERRNM(merr);
1212: #endif
1213: }
1214: PetscCall(PetscInfo(NULL, "Found %zu boundary vertices, %zu boundary faces and %zu boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size()));
1216: /* Get the material sets and populate the data for all locally owned elements */
1217: {
1218: PetscCall(PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials));
1219: /* Get the count of entities of particular type from dmmoab->elocal
1220: -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1221: moab::Range msets;
1222: merr = dmmoab->mbiface->get_entities_by_type_and_tag(dmmoab->fileset, moab::MBENTITYSET, &dmmoab->material_tag, NULL, 1, msets, moab::Interface::UNION);
1223: MBERRNM(merr);
1224: if (msets.size() == 0) PetscCall(PetscInfo(NULL, "No material sets found in the fileset.\n"));
1226: for (unsigned i = 0; i < msets.size(); ++i) {
1227: moab::Range msetelems;
1228: merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);
1229: MBERRNM(merr);
1230: #ifdef MOAB_HAVE_MPI
1231: /* filter all the non-owned and shared entities out of the list */
1232: merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1233: MBERRNM(merr);
1234: #endif
1236: int partID;
1237: moab::EntityHandle mset = msets[i];
1238: merr = dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID);
1239: MBERRNM(merr);
1241: for (unsigned j = 0; j < msetelems.size(); ++j) dmmoab->materials[dmmoab->elocal->index(msetelems[j])] = partID;
1242: }
1243: }
1244: PetscFunctionReturn(PETSC_SUCCESS);
1245: }
1247: /*@C
1248: DMMoabCreateVertices - Creates and adds several vertices to the primary set represented by the DM.
1250: Collective
1252: Input Parameters:
1253: + dm - The `DM` object
1254: . coords - The connectivity of the element
1255: - nverts - The number of vertices that form the element
1257: Output Parameter:
1258: . overts - The list of vertices that were created (can be `NULL`)
1260: Level: beginner
1262: .seealso: `DMMOAB`, `DMMoabCreateSubmesh()`, `DMMoabCreateElement()`
1263: @*/
1264: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal *coords, PetscInt nverts, moab::Range *overts)
1265: {
1266: moab::ErrorCode merr;
1267: DM_Moab *dmmoab;
1268: moab::Range verts;
1270: PetscFunctionBegin;
1272: PetscAssertPointer(coords, 2);
1274: dmmoab = (DM_Moab *)dm->data;
1276: /* Insert new points */
1277: merr = dmmoab->mbiface->create_vertices(&coords[0], nverts, verts);
1278: MBERRNM(merr);
1279: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, verts);
1280: MBERRNM(merr);
1282: if (overts) *overts = verts;
1283: PetscFunctionReturn(PETSC_SUCCESS);
1284: }
1286: /*@C
1287: DMMoabCreateElement - Adds an element of specified type to the primary set represented by the DM.
1289: Collective
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: . oelem - The handle to the element created and added to the `DM` object
1300: Level: beginner
1302: .seealso: `DMMOAB`, `DMMoabCreateSubmesh()`, `DMMoabCreateVertices()`
1303: @*/
1304: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle *conn, PetscInt nverts, moab::EntityHandle *oelem)
1305: {
1306: moab::ErrorCode merr;
1307: DM_Moab *dmmoab;
1308: moab::EntityHandle elem;
1310: PetscFunctionBegin;
1312: PetscAssertPointer(conn, 3);
1314: dmmoab = (DM_Moab *)dm->data;
1316: /* Insert new element */
1317: merr = dmmoab->mbiface->create_element(type, conn, nverts, elem);
1318: MBERRNM(merr);
1319: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1);
1320: MBERRNM(merr);
1322: if (oelem) *oelem = elem;
1323: PetscFunctionReturn(PETSC_SUCCESS);
1324: }
1326: /*@C
1327: DMMoabCreateSubmesh - Creates a sub-`DM` object with a set that contains all vertices/elements of the parent
1328: in addition to providing support for dynamic mesh modifications. This is useful for AMR calculations to
1329: create a DM object on a refined level.
1331: Collective
1333: Input Parameters:
1334: . dm - The `DM` object
1336: Output Parameter:
1337: . newdm - The sub `DM` object with updated set information
1339: Level: advanced
1341: .seealso: `DMMOAB`, `DMCreate()`, `DMMoabCreateVertices()`, `DMMoabCreateElement()`
1342: @*/
1343: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1344: {
1345: DM_Moab *dmmoab;
1346: DM_Moab *ndmmoab;
1347: moab::ErrorCode merr;
1349: PetscFunctionBegin;
1352: dmmoab = (DM_Moab *)dm->data;
1354: /* Create the basic DMMOAB object and keep the default parameters created by DM impls */
1355: PetscCall(DMMoabCreateMoab(((PetscObject)dm)->comm, dmmoab->mbiface, &dmmoab->ltog_tag, NULL, newdm));
1357: /* get all the necessary handles from the private DM object */
1358: ndmmoab = (DM_Moab *)(*newdm)->data;
1360: /* set the sub-mesh's parent DM reference */
1361: ndmmoab->parent = &dm;
1363: /* create a file set to associate all entities in current mesh */
1364: merr = ndmmoab->mbiface->create_meshset(moab::MESHSET_SET, ndmmoab->fileset);
1365: MBERR("Creating file set failed", merr);
1367: /* create a meshset and then add old fileset as child */
1368: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal);
1369: MBERR("Adding child vertices to parent failed", merr);
1370: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal);
1371: MBERR("Adding child elements to parent failed", merr);
1373: /* preserve the field association between the parent and sub-mesh objects */
1374: PetscCall(DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames));
1375: PetscFunctionReturn(PETSC_SUCCESS);
1376: }
1378: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1379: {
1380: DM_Moab *dmmoab = (DM_Moab *)dm->data;
1381: const char *name;
1382: MPI_Comm comm;
1383: PetscMPIInt size;
1385: PetscFunctionBegin;
1386: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1387: PetscCallMPI(MPI_Comm_size(comm, &size));
1388: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1389: PetscCall(PetscViewerASCIIPushTab(viewer));
1390: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimensions:\n", name, dmmoab->dim));
1391: else PetscCall(PetscViewerASCIIPrintf(viewer, "Mesh in %" PetscInt_FMT " dimensions:\n", dmmoab->dim));
1392: /* print details about the global mesh */
1393: {
1394: PetscCall(PetscViewerASCIIPushTab(viewer));
1395: PetscCall(PetscViewerASCIIPrintf(viewer, "Sizes: cells=%" PetscInt_FMT ", vertices=%" PetscInt_FMT ", blocks=%" PetscInt_FMT "\n", dmmoab->nele, dmmoab->n, dmmoab->bs));
1396: /* print boundary data */
1397: PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary trace:\n"));
1398: {
1399: PetscCall(PetscViewerASCIIPushTab(viewer));
1400: PetscCall(PetscViewerASCIIPrintf(viewer, "cells=%zu, faces=%zu, vertices=%zu\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size()));
1401: PetscCall(PetscViewerASCIIPopTab(viewer));
1402: }
1403: /* print field data */
1404: PetscCall(PetscViewerASCIIPrintf(viewer, "Fields: %" PetscInt_FMT " components\n", dmmoab->numFields));
1405: {
1406: PetscCall(PetscViewerASCIIPushTab(viewer));
1407: for (int i = 0; i < dmmoab->numFields; ++i) PetscCall(PetscViewerASCIIPrintf(viewer, "[%" PetscInt_FMT "] - %s\n", i, dmmoab->fieldNames[i]));
1408: PetscCall(PetscViewerASCIIPopTab(viewer));
1409: }
1410: PetscCall(PetscViewerASCIIPopTab(viewer));
1411: }
1412: PetscCall(PetscViewerASCIIPopTab(viewer));
1413: PetscCall(PetscViewerFlush(viewer));
1414: PetscFunctionReturn(PETSC_SUCCESS);
1415: }
1417: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1418: {
1419: PetscFunctionReturn(PETSC_SUCCESS);
1420: }
1422: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1423: {
1424: PetscFunctionReturn(PETSC_SUCCESS);
1425: }
1427: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1428: {
1429: PetscBool isascii, ishdf5, isvtk;
1431: PetscFunctionBegin;
1434: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1435: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1436: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1437: if (isascii) {
1438: PetscCall(DMMoabView_Ascii(dm, viewer));
1439: } else if (ishdf5) {
1440: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1441: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ));
1442: PetscCall(DMMoabView_HDF5(dm, viewer));
1443: PetscCall(PetscViewerPopFormat(viewer));
1444: #else
1445: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1446: #endif
1447: } else if (isvtk) PetscCall(DMMoabView_VTK(dm, viewer));
1448: PetscFunctionReturn(PETSC_SUCCESS);
1449: }
1451: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
1452: {
1453: PetscFunctionBegin;
1454: dm->ops->view = DMView_Moab;
1455: dm->ops->load = NULL /* DMLoad_Moab */;
1456: dm->ops->setfromoptions = DMSetFromOptions_Moab;
1457: dm->ops->clone = DMClone_Moab;
1458: dm->ops->setup = DMSetUp_Moab;
1459: dm->ops->createlocalsection = NULL;
1460: dm->ops->createsectionpermutation = NULL;
1461: dm->ops->createdefaultconstraints = NULL;
1462: dm->ops->createglobalvector = DMCreateGlobalVector_Moab;
1463: dm->ops->createlocalvector = DMCreateLocalVector_Moab;
1464: dm->ops->getlocaltoglobalmapping = NULL;
1465: dm->ops->createfieldis = NULL;
1466: dm->ops->createcoordinatedm = NULL /* DMCreateCoordinateDM_Moab */;
1467: dm->ops->createcellcoordinatedm = NULL;
1468: dm->ops->getcoloring = NULL;
1469: dm->ops->creatematrix = DMCreateMatrix_Moab;
1470: dm->ops->createinterpolation = DMCreateInterpolation_Moab;
1471: dm->ops->createinjection = NULL /* DMCreateInjection_Moab */;
1472: dm->ops->refine = DMRefine_Moab;
1473: dm->ops->coarsen = DMCoarsen_Moab;
1474: dm->ops->refinehierarchy = DMRefineHierarchy_Moab;
1475: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Moab;
1476: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Moab;
1477: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Moab;
1478: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Moab;
1479: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Moab;
1480: dm->ops->destroy = DMDestroy_Moab;
1481: dm->ops->createsubdm = NULL /* DMCreateSubDM_Moab */;
1482: dm->ops->getdimpoints = NULL /* DMGetDimPoints_Moab */;
1483: dm->ops->locatepoints = NULL /* DMLocatePoints_Moab */;
1484: PetscFunctionReturn(PETSC_SUCCESS);
1485: }
1487: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1488: {
1489: PetscFunctionBegin;
1490: /* get all the necessary handles from the private DM object */
1491: (*newdm)->data = (DM_Moab *)dm->data;
1492: ((DM_Moab *)dm->data)->refct++;
1494: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMMOAB));
1495: PetscCall(DMInitialize_Moab(*newdm));
1496: PetscFunctionReturn(PETSC_SUCCESS);
1497: }
1499: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1500: {
1501: PetscFunctionBegin;
1503: PetscCall(PetscNew((DM_Moab **)&dm->data));
1505: ((DM_Moab *)dm->data)->bs = 1;
1506: ((DM_Moab *)dm->data)->numFields = 1;
1507: ((DM_Moab *)dm->data)->n = 0;
1508: ((DM_Moab *)dm->data)->nloc = 0;
1509: ((DM_Moab *)dm->data)->nghost = 0;
1510: ((DM_Moab *)dm->data)->nele = 0;
1511: ((DM_Moab *)dm->data)->neleloc = 0;
1512: ((DM_Moab *)dm->data)->neleghost = 0;
1513: ((DM_Moab *)dm->data)->ltog_map = NULL;
1514: ((DM_Moab *)dm->data)->ltog_sendrecv = NULL;
1516: ((DM_Moab *)dm->data)->refct = 1;
1517: ((DM_Moab *)dm->data)->parent = NULL;
1518: ((DM_Moab *)dm->data)->vlocal = new moab::Range();
1519: ((DM_Moab *)dm->data)->vowned = new moab::Range();
1520: ((DM_Moab *)dm->data)->vghost = new moab::Range();
1521: ((DM_Moab *)dm->data)->elocal = new moab::Range();
1522: ((DM_Moab *)dm->data)->eghost = new moab::Range();
1524: PetscCall(DMInitialize_Moab(dm));
1525: PetscFunctionReturn(PETSC_SUCCESS);
1526: }