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_INTERN PetscErrorCode DMCreateInterpolation_Moab(DM, DM, Mat *, Vec *);
24: PETSC_INTERN PetscErrorCode DMCreateMatrix_Moab(DM, Mat *);
25: PETSC_INTERN PetscErrorCode DMRefine_Moab(DM, MPI_Comm, DM *);
26: PETSC_INTERN PetscErrorCode DMCoarsen_Moab(DM, MPI_Comm, DM *);
27: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Moab(DM, PetscInt, DM[]);
28: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM, PetscInt, DM[]);
29: PETSC_INTERN PetscErrorCode DMCreateGlobalVector_Moab(DM, Vec *);
30: PETSC_INTERN PetscErrorCode DMCreateLocalVector_Moab(DM, Vec *);
31: PETSC_INTERN PetscErrorCode DMCreateMatrix_Moab(DM, Mat *J);
32: PETSC_INTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM, Vec, InsertMode, Vec);
33: PETSC_INTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM, Vec, InsertMode, Vec);
34: PETSC_INTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM, Vec, InsertMode, Vec);
35: PETSC_INTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM, Vec, InsertMode, Vec);
37: /* Un-implemented routines */
38: /*
39: PETSC_INTERN PetscErrorCode DMCreatelocalsection_Moab(DM);
40: PETSC_INTERN PetscErrorCode DMCreateInjection_Moab(DM, DM, Mat *);
41: PETSC_INTERN PetscErrorCode DMLoad_Moab(DM, PetscViewer);
42: PETSC_INTERN PetscErrorCode DMGetDimPoints_Moab(DM, PetscInt, PetscInt *, PetscInt *);
43: PETSC_INTERN PetscErrorCode DMCreateSubDM_Moab(DM, PetscInt, PetscInt [], IS *, DM *);
44: PETSC_INTERN PetscErrorCode DMLocatePoints_Moab(DM, Vec, IS *);
45: */
47: /*@C
48: DMMoabCreate - Creates a `DMMOAB` object, which encapsulates a moab instance
50: Collective
52: Input Parameter:
53: . comm - The communicator for the `DMMOAB` object
55: Output Parameter:
56: . dmb - The `DMMOAB` object
58: Level: beginner
60: .seealso: `DMMOAB`, `DMMoabCreateMoab()`
61: @*/
62: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
63: {
64: PetscFunctionBegin;
65: PetscAssertPointer(dmb, 2);
66: PetscCall(DMCreate(comm, dmb));
67: PetscCall(DMSetType(*dmb, DMMOAB));
68: PetscFunctionReturn(PETSC_SUCCESS);
69: }
71: /*@C
72: DMMoabCreateMoab - Creates a `DMMOAB` object, optionally from an instance and other data
74: Collective
76: Input Parameters:
77: + comm - The communicator for the `DMMOAB` object
78: . mbiface - (ptr to) the `DMMOAB` Instance; if passed in `NULL`, MOAB instance is created inside PETSc, and destroyed
79: along with the `DMMOAB`
80: . ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
81: - range - If non-`NULL`, contains range of entities to which DOFs will be assigned
83: Output Parameter:
84: . dmb - The `DMMOAB` object
86: Level: intermediate
88: .seealso: `DMMOAB`, `DMMoabCreate()`
89: @*/
90: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
91: {
92: DM dmmb;
93: DM_Moab *dmmoab;
95: PetscFunctionBegin;
96: PetscAssertPointer(dmb, 6);
98: PetscCall(DMMoabCreate(comm, &dmmb));
99: dmmoab = (DM_Moab *)dmmb->data;
101: if (!mbiface) {
102: dmmoab->mbiface = new moab::Core();
103: dmmoab->icreatedinstance = PETSC_TRUE;
104: } else {
105: dmmoab->mbiface = mbiface;
106: dmmoab->icreatedinstance = PETSC_FALSE;
107: }
109: /* by default the fileset = root set. This set stores the hierarchy of entities belonging to current DM */
110: dmmoab->fileset = 0;
111: dmmoab->hlevel = 0;
112: dmmoab->nghostrings = 0;
114: #if defined(MOAB_HAVE_MPI)
115: moab::EntityHandle partnset;
117: /* Create root sets for each mesh. Then pass these
118: to the load_file functions to be populated. */
119: PetscCallMOAB(dmmoab->mbiface->create_meshset(moab::MESHSET_SET, partnset));
121: /* Create the parallel communicator object with the partition handle associated with MOAB */
122: dmmoab->pcomm = moab::ParallelComm::get_pcomm(dmmoab->mbiface, partnset, &comm);
123: #endif
125: /* do the remaining initializations for DMMoab */
126: dmmoab->bs = 1;
127: dmmoab->numFields = 1;
128: PetscCall(PetscMalloc(dmmoab->numFields * sizeof(char *), &dmmoab->fieldNames));
129: PetscCall(PetscStrallocpy("DEFAULT", (char **)&dmmoab->fieldNames[0]));
130: dmmoab->rw_dbglevel = 0;
131: dmmoab->partition_by_rank = PETSC_FALSE;
132: dmmoab->extra_read_options[0] = '\0';
133: dmmoab->extra_write_options[0] = '\0';
134: dmmoab->read_mode = READ_PART;
135: dmmoab->write_mode = WRITE_PART;
137: /* set global ID tag handle */
138: if (ltog_tag && *ltog_tag) PetscCall(DMMoabSetLocalToGlobalTag(dmmb, *ltog_tag));
139: else {
140: PetscCallMOAB(dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag));
141: if (ltog_tag) *ltog_tag = dmmoab->ltog_tag;
142: }
144: PetscCallMOAB(dmmoab->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dmmoab->material_tag));
146: /* set the local range of entities (vertices) of interest */
147: if (range) PetscCall(DMMoabSetLocalVertices(dmmb, range));
148: *dmb = dmmb;
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: #if defined(MOAB_HAVE_MPI)
154: /*@C
155: DMMoabGetParallelComm - Get the ParallelComm used with this `DMMOAB`
157: Collective
159: Input Parameter:
160: . dm - The `DMMOAB` object being set
162: Output Parameter:
163: . pcomm - The ParallelComm for the `DMMOAB`
165: Level: beginner
167: .seealso: `DMMOAB`, `DMMoabSetInterface()`
168: @*/
169: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
170: {
171: PetscFunctionBegin;
173: *pcomm = ((DM_Moab *)dm->data)->pcomm;
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: #endif /* MOAB_HAVE_MPI */
179: /*@C
180: DMMoabSetInterface - Set the MOAB instance used with this `DMMOAB`
182: Collective
184: Input Parameters:
185: + dm - The `DMMOAB` object being set
186: - mbiface - The MOAB instance being set on this `DMMOAB`
188: Level: beginner
190: .seealso: `DMMOAB`, `DMMoabGetInterface()`
191: @*/
192: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
193: {
194: DM_Moab *dmmoab = (DM_Moab *)dm->data;
196: PetscFunctionBegin;
198: PetscAssertPointer(mbiface, 2);
199: #if defined(MOAB_HAVE_MPI)
200: dmmoab->pcomm = NULL;
201: #endif
202: dmmoab->mbiface = mbiface;
203: dmmoab->icreatedinstance = PETSC_FALSE;
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: /*@C
208: DMMoabGetInterface - Get the MOAB instance used with this `DMMOAB`
210: Collective
212: Input Parameter:
213: . dm - The `DMMOAB` object being set
215: Output Parameter:
216: . mbiface - The MOAB instance set on this `DMMOAB`
218: Level: beginner
220: .seealso: `DMMOAB`, `DMMoabSetInterface()`
221: @*/
222: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
223: {
224: static PetscBool cite = PETSC_FALSE;
226: PetscFunctionBegin;
228: PetscCall(PetscCitationsRegister("@techreport{tautges_moab:_2004,\n type = {{SAND2004-1592}},\n title = {{MOAB:} A Mesh-Oriented Database}, institution = {Sandia National Laboratories},\n author = {Tautges, T. J. and Meyers, R. and Merkley, "
229: "K. and Stimpson, C. and Ernst, C.},\n year = {2004}, note = {Report}\n}\n",
230: &cite));
231: *mbiface = ((DM_Moab *)dm->data)->mbiface;
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: /*@C
236: DMMoabSetLocalVertices - Set the entities having DOFs on this `DMMOAB`
238: Collective
240: Input Parameters:
241: + dm - The `DMMOAB` object being set
242: - range - The entities treated by this `DMMOAB`
244: Level: beginner
246: .seealso: `DMMOAB`, `DMMoabGetAllVertices()`
247: @*/
248: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
249: {
250: moab::Range tmpvtxs;
251: DM_Moab *dmmoab = (DM_Moab *)dm->data;
253: PetscFunctionBegin;
255: dmmoab->vlocal->clear();
256: dmmoab->vowned->clear();
258: dmmoab->vlocal->insert(range->begin(), range->end());
260: #if defined(MOAB_HAVE_MPI)
261: /* filter based on parallel status */
262: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned));
264: /* filter all the non-owned and shared entities out of the list */
265: tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
266: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost));
267: tmpvtxs = moab::subtract(tmpvtxs, *dmmoab->vghost);
268: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
269: #else
270: *dmmoab->vowned = *dmmoab->vlocal;
271: #endif
273: /* compute and cache the sizes of local and ghosted entities */
274: dmmoab->nloc = dmmoab->vowned->size();
275: dmmoab->nghost = dmmoab->vghost->size();
276: #if defined(MOAB_HAVE_MPI)
277: PetscCallMPI(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
278: #else
279: dmmoab->n = dmmoab->nloc;
280: #endif
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }
284: /*@C
285: DMMoabGetAllVertices - Get the entities having DOFs on this `DMMOAB`
287: Collective
289: Input Parameter:
290: . dm - The `DMMOAB` object being set
292: Output Parameter:
293: . local - The local vertex entities in this `DMMOAB` = (owned+ghosted)
295: Level: beginner
297: .seealso: `DMMOAB`, `DMMoabGetLocalVertices()`
298: @*/
299: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
300: {
301: PetscFunctionBegin;
303: if (local) *local = *((DM_Moab *)dm->data)->vlocal;
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }
307: /*@C
308: DMMoabGetLocalVertices - Get the entities having DOFs on this `DMMOAB`
310: Collective
312: Input Parameter:
313: . dm - The `DMMOAB` object being set
315: Output Parameters:
316: + owned - The owned vertex entities in this `DMMOAB`
317: - ghost - The ghosted entities (non-owned) stored locally in this partition
319: Level: beginner
321: .seealso: `DMMOAB`, `DMMoabGetAllVertices()`
322: @*/
323: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
324: {
325: PetscFunctionBegin;
327: if (owned) *owned = ((DM_Moab *)dm->data)->vowned;
328: if (ghost) *ghost = ((DM_Moab *)dm->data)->vghost;
329: PetscFunctionReturn(PETSC_SUCCESS);
330: }
332: /*@C
333: DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned
335: Collective
337: Input Parameter:
338: . dm - The `DMMOAB` object being set
340: Output Parameter:
341: . range - The entities owned locally
343: Level: beginner
345: .seealso: `DMMOAB`, `DMMoabSetLocalElements()`
346: @*/
347: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
348: {
349: PetscFunctionBegin;
351: if (range) *range = ((DM_Moab *)dm->data)->elocal;
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /*@C
356: DMMoabSetLocalElements - Set the entities having DOFs on this `DMMOAB`
358: Collective
360: Input Parameters:
361: + dm - The `DMMOAB` object being set
362: - range - The entities treated by this `DMMOAB`
364: Level: beginner
366: .seealso: `DMMOAB`, `DMMoabGetLocalElements()`
367: @*/
368: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
369: {
370: DM_Moab *dmmoab = (DM_Moab *)dm->data;
372: PetscFunctionBegin;
374: dmmoab->elocal->clear();
375: dmmoab->eghost->clear();
376: dmmoab->elocal->insert(range->begin(), range->end());
377: #if defined(MOAB_HAVE_MPI)
378: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT));
379: *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
380: #endif
381: dmmoab->neleloc = dmmoab->elocal->size();
382: dmmoab->neleghost = dmmoab->eghost->size();
383: #if defined(MOAB_HAVE_MPI)
384: PetscCallMPI(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
385: PetscCall(PetscInfo(dm, "Created %" PetscInt_FMT " local and %" PetscInt_FMT " global elements.\n", dmmoab->neleloc, dmmoab->nele));
386: #else
387: dmmoab->nele = dmmoab->neleloc;
388: #endif
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: /*@C
393: DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
395: Collective
397: Input Parameters:
398: + dm - The `DMMOAB` object being set
399: - ltogtag - The `DMMOAB` tag used for local to global ids
401: Level: beginner
403: .seealso: `DMMOAB`, `DMMoabGetLocalToGlobalTag()`
404: @*/
405: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
406: {
407: PetscFunctionBegin;
409: ((DM_Moab *)dm->data)->ltog_tag = ltogtag;
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: /*@C
414: DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
416: Collective
418: Input Parameter:
419: . dm - The `DMMOAB` object being set
421: Output Parameter:
422: . ltog_tag - The MOAB tag used for local to global ids
424: Level: beginner
426: .seealso: `DMMOAB`, `DMMoabSetLocalToGlobalTag()`
427: @*/
428: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
429: {
430: PetscFunctionBegin;
432: *ltog_tag = ((DM_Moab *)dm->data)->ltog_tag;
433: PetscFunctionReturn(PETSC_SUCCESS);
434: }
436: /*@C
437: DMMoabSetBlockSize - Set the block size used with this `DMMOAB`
439: Collective
441: Input Parameters:
442: + dm - The `DMMOAB` object being set
443: - bs - The block size used with this `DMMOAB`
445: Level: beginner
447: .seealso: `DMMOAB`, `DMMoabGetBlockSize()`
448: @*/
449: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
450: {
451: PetscFunctionBegin;
453: ((DM_Moab *)dm->data)->bs = bs;
454: PetscFunctionReturn(PETSC_SUCCESS);
455: }
457: /*@C
458: DMMoabGetBlockSize - Get the block size used with this `DMMOAB`
460: Collective
462: Input Parameter:
463: . dm - The `DMMOAB` object being set
465: Output Parameter:
466: . bs - The block size used with this `DMMOAB`
468: Level: beginner
470: .seealso: `DMMOAB`, `DMMoabSetBlockSize()`
471: @*/
472: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
473: {
474: PetscFunctionBegin;
476: *bs = ((DM_Moab *)dm->data)->bs;
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }
480: /*@C
481: DMMoabGetSize - Get the global vertex size used with this `DMMOAB`
483: Collective
485: Input Parameter:
486: . dm - The `DMMOAB` object being set
488: Output Parameters:
489: + neg - The number of global elements in the `DMMOAB` instance
490: - nvg - The number of global vertices in the `DMMOAB` instance
492: Level: beginner
494: .seealso: `DMMOAB`, `DMMoabGetLocalSize()`
495: @*/
496: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
497: {
498: PetscFunctionBegin;
500: if (neg) *neg = ((DM_Moab *)dm->data)->nele;
501: if (nvg) *nvg = ((DM_Moab *)dm->data)->n;
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }
505: /*@C
506: DMMoabGetLocalSize - Get the local and ghosted vertex size used with this `DMMOAB`
508: Collective
510: Input Parameter:
511: . dm - The `DMMOAB` object being set
513: Output Parameters:
514: + nel - The number of owned elements in this processor
515: . neg - The number of ghosted elements in this processor
516: . nvl - The number of owned vertices in this processor
517: - nvg - The number of ghosted vertices in this processor
519: Level: beginner
521: .seealso: `DMMOAB`, `DMMoabGetSize()`
522: @*/
523: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
524: {
525: PetscFunctionBegin;
527: if (nel) *nel = ((DM_Moab *)dm->data)->neleloc;
528: if (neg) *neg = ((DM_Moab *)dm->data)->neleghost;
529: if (nvl) *nvl = ((DM_Moab *)dm->data)->nloc;
530: if (nvg) *nvg = ((DM_Moab *)dm->data)->nghost;
531: PetscFunctionReturn(PETSC_SUCCESS);
532: }
534: /*@C
535: DMMoabGetOffset - Get the local offset for the global vector
537: Collective
539: Input Parameter:
540: . dm - The `DMMOAB` object being set
542: Output Parameter:
543: . offset - The local offset for the global vector
545: Level: beginner
547: .seealso: `DMMOAB`, `DMMoabGetDimension()`
548: @*/
549: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
550: {
551: PetscFunctionBegin;
553: *offset = ((DM_Moab *)dm->data)->vstart;
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }
557: /*@C
558: DMMoabGetDimension - Get the dimension of the `DM` Mesh
560: Collective
562: Input Parameter:
563: . dm - The `DMMOAB` object
565: Output Parameter:
566: . dim - The dimension of `DM`
568: Level: beginner
570: .seealso: `DMMOAB`, `DMMoabGetOffset()`
571: @*/
572: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
573: {
574: PetscFunctionBegin;
576: *dim = ((DM_Moab *)dm->data)->dim;
577: PetscFunctionReturn(PETSC_SUCCESS);
578: }
580: /*@C
581: DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
582: generated through uniform refinement.
584: Collective
586: Input Parameter:
587: . dm - The `DMMOAB` object being set
589: Output Parameter:
590: . nlevel - The current mesh hierarchy level
592: Level: beginner
594: .seealso: `DMMOAB`, `DMMoabGetMaterialBlock()`
595: @*/
596: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
597: {
598: PetscFunctionBegin;
600: if (nlevel) *nlevel = ((DM_Moab *)dm->data)->hlevel;
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
604: /*@C
605: DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh
607: Collective
609: Input Parameters:
610: + dm - The `DMMOAB` object
611: - ehandle - The element entity handle
613: Output Parameter:
614: . mat - The material ID for the current entity
616: Level: beginner
618: .seealso: `DMMOAB`, `DMMoabGetHierarchyLevel()`
619: @*/
620: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
621: {
622: DM_Moab *dmmoab;
624: PetscFunctionBegin;
626: if (*mat) {
627: dmmoab = (DM_Moab *)dm->data;
628: *mat = dmmoab->materials[dmmoab->elocal->index(ehandle)];
629: }
630: PetscFunctionReturn(PETSC_SUCCESS);
631: }
633: /*@C
634: DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities
636: Collective
638: Input Parameters:
639: + dm - The `DMMOAB` object
640: . nconn - Number of entities whose coordinates are needed
641: - conn - The vertex entity handles
643: Output Parameter:
644: . vpos - The coordinates of the requested vertex entities
646: Level: beginner
648: .seealso: `DMMOAB`, `DMMoabGetVertexConnectivity()`
649: @*/
650: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
651: {
652: DM_Moab *dmmoab;
654: PetscFunctionBegin;
656: PetscAssertPointer(conn, 3);
657: PetscAssertPointer(vpos, 4);
658: dmmoab = (DM_Moab *)dm->data;
660: /* Get connectivity information in MOAB canonical ordering */
661: if (dmmoab->hlevel) PetscCallMOAB(dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle *>(conn), nconn, dmmoab->hlevel, vpos));
662: else PetscCallMOAB(dmmoab->mbiface->get_coords(conn, nconn, vpos));
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: /*@C
667: DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity
669: Collective
671: Input Parameters:
672: + dm - The `DMMOAB` object
673: - vhandle - Vertex entity handle
675: Output Parameters:
676: + nconn - Number of entities whose coordinates are needed
677: - conn - The vertex entity handles
679: Level: beginner
681: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabRestoreVertexConnectivity()`
682: @*/
683: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt *nconn, moab::EntityHandle **conn)
684: {
685: DM_Moab *dmmoab;
686: std::vector<moab::EntityHandle> adj_entities, connect;
688: PetscFunctionBegin;
690: PetscAssertPointer(conn, 4);
691: dmmoab = (DM_Moab *)dm->data;
693: /* Get connectivity information in MOAB canonical ordering */
694: PetscCallMOAB(dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION));
695: PetscCallMOAB(dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect));
697: if (conn) {
698: PetscCall(PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn));
699: PetscCall(PetscArraycpy(*conn, &connect[0], connect.size()));
700: }
701: if (nconn) *nconn = connect.size();
702: PetscFunctionReturn(PETSC_SUCCESS);
703: }
705: /*@C
706: DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity
708: Collective
710: Input Parameters:
711: + dm - The `DMMOAB` object
712: . ehandle - Vertex entity handle
713: . nconn - Number of entities whose coordinates are needed
714: - conn - The vertex entity handles
716: Level: beginner
718: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`
719: @*/
720: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, moab::EntityHandle **conn)
721: {
722: PetscFunctionBegin;
724: PetscAssertPointer(conn, 4);
726: if (conn) PetscCall(PetscFree(*conn));
727: if (nconn) *nconn = 0;
728: PetscFunctionReturn(PETSC_SUCCESS);
729: }
731: /*@C
732: DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity
734: Collective
736: Input Parameters:
737: + dm - The `DMMOAB` object
738: - ehandle - Vertex entity handle
740: Output Parameters:
741: + nconn - Number of entities whose coordinates are needed
742: - conn - The vertex entity handles
744: Level: beginner
746: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`, `DMMoabRestoreVertexConnectivity()`
747: @*/
748: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, const moab::EntityHandle **conn)
749: {
750: DM_Moab *dmmoab;
751: const moab::EntityHandle *connect;
752: std::vector<moab::EntityHandle> vconn;
753: PetscInt nnodes;
755: PetscFunctionBegin;
757: PetscAssertPointer(conn, 4);
758: dmmoab = (DM_Moab *)dm->data;
760: /* Get connectivity information in MOAB canonical ordering */
761: PetscCallMOAB(dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes));
762: if (conn) *conn = connect;
763: if (nconn) *nconn = nnodes;
764: PetscFunctionReturn(PETSC_SUCCESS);
765: }
767: /*@C
768: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
770: Collective
772: Input Parameters:
773: + dm - The `DMMOAB` object
774: - ent - Entity handle
776: Output Parameter:
777: . ent_on_boundary - `PETSC_TRUE` if entity on boundary; `PETSC_FALSE` otherwise
779: Level: beginner
781: .seealso: `DMMOAB`, `DMMoabCheckBoundaryVertices()`
782: @*/
783: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool *ent_on_boundary)
784: {
785: moab::EntityType etype;
786: DM_Moab *dmmoab;
787: PetscInt edim;
789: PetscFunctionBegin;
791: PetscAssertPointer(ent_on_boundary, 3);
792: dmmoab = (DM_Moab *)dm->data;
794: /* get the entity type and handle accordingly */
795: etype = dmmoab->mbiface->type_from_handle(ent);
796: PetscCheck(etype < moab::MBPOLYHEDRON, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Entity type on the boundary skin is invalid. EntityType = %" PetscInt_FMT, etype);
798: /* get the entity dimension */
799: edim = dmmoab->mbiface->dimension_from_handle(ent);
801: *ent_on_boundary = PETSC_FALSE;
802: if (etype == moab::MBVERTEX && edim == 0) {
803: *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
804: } else {
805: if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
806: if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
807: } else { /* next check the lower-dimensional faces */
808: if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
809: }
810: }
811: PetscFunctionReturn(PETSC_SUCCESS);
812: }
814: /*@C
815: DMMoabCheckBoundaryVertices - Check whether a given entity is on the boundary (vertex, edge, face, element)
817: Input Parameters:
818: + dm - The `DMMOAB` object
819: . nconn - Number of handles
820: - cnt - Array of entity handles
822: Output Parameter:
823: . isbdvtx - Array of boundary markers - `PETSC_TRUE` if entity on boundary; `PETSC_FALSE` otherwise
825: Level: beginner
827: .seealso: `DMMOAB`, `DMMoabIsEntityOnBoundary()`
828: @*/
829: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool *isbdvtx)
830: {
831: DM_Moab *dmmoab;
832: PetscInt i;
834: PetscFunctionBegin;
836: PetscAssertPointer(cnt, 3);
837: PetscAssertPointer(isbdvtx, 4);
838: dmmoab = (DM_Moab *)dm->data;
840: for (i = 0; i < nconn; ++i) isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
841: PetscFunctionReturn(PETSC_SUCCESS);
842: }
844: /*@C
845: DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary
847: Input Parameter:
848: . dm - The `DMMOAB` object
850: Output Parameters:
851: + bdvtx - Boundary vertices
852: . bdelems - Boundary elements
853: - bdfaces - Boundary faces
855: Level: beginner
857: .seealso: `DMMOAB`, `DMMoabCheckBoundaryVertices()`, `DMMoabIsEntityOnBoundary()`
858: @*/
859: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range **bdelems, const moab::Range **bdfaces)
860: {
861: DM_Moab *dmmoab;
863: PetscFunctionBegin;
865: dmmoab = (DM_Moab *)dm->data;
867: if (bdvtx) *bdvtx = dmmoab->bndyvtx;
868: if (bdfaces) *bdfaces = dmmoab->bndyfaces;
869: if (bdelems) *bdfaces = dmmoab->bndyelems;
870: PetscFunctionReturn(PETSC_SUCCESS);
871: }
873: static PetscErrorCode DMDestroy_Moab(DM dm)
874: {
875: PetscInt i;
876: DM_Moab *dmmoab = (DM_Moab *)dm->data;
878: PetscFunctionBegin;
881: dmmoab->refct--;
882: if (!dmmoab->refct) {
883: delete dmmoab->vlocal;
884: delete dmmoab->vowned;
885: delete dmmoab->vghost;
886: delete dmmoab->elocal;
887: delete dmmoab->eghost;
888: delete dmmoab->bndyvtx;
889: delete dmmoab->bndyfaces;
890: delete dmmoab->bndyelems;
892: PetscCall(PetscFree(dmmoab->gsindices));
893: PetscCall(PetscFree2(dmmoab->gidmap, dmmoab->lidmap));
894: PetscCall(PetscFree(dmmoab->dfill));
895: PetscCall(PetscFree(dmmoab->ofill));
896: PetscCall(PetscFree(dmmoab->materials));
897: if (dmmoab->fieldNames) {
898: for (i = 0; i < dmmoab->numFields; i++) PetscCall(PetscFree(dmmoab->fieldNames[i]));
899: PetscCall(PetscFree(dmmoab->fieldNames));
900: }
902: if (dmmoab->nhlevels) {
903: PetscCall(PetscFree(dmmoab->hsets));
904: dmmoab->nhlevels = 0;
905: if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
906: dmmoab->hierarchy = NULL;
907: }
909: if (dmmoab->icreatedinstance) {
910: delete dmmoab->pcomm;
911: PetscCallMOAB(dmmoab->mbiface->delete_mesh());
912: delete dmmoab->mbiface;
913: }
914: dmmoab->mbiface = NULL;
915: #if defined(MOAB_HAVE_MPI)
916: dmmoab->pcomm = NULL;
917: #endif
918: PetscCall(VecScatterDestroy(&dmmoab->ltog_sendrecv));
919: PetscCall(ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map));
920: PetscCall(PetscFree(dm->data));
921: }
922: PetscFunctionReturn(PETSC_SUCCESS);
923: }
925: static PetscErrorCode DMSetFromOptions_Moab(DM dm, PetscOptionItems PetscOptionsObject)
926: {
927: DM_Moab *dmmoab = (DM_Moab *)dm->data;
929: PetscFunctionBegin;
930: PetscOptionsHeadBegin(PetscOptionsObject, "DMMoab Options");
931: PetscCall(PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL, 0));
932: 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));
933: /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
934: 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));
935: 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));
936: PetscCall(PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum *)&dmmoab->read_mode, NULL));
937: PetscCall(PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum *)&dmmoab->write_mode, NULL));
938: PetscOptionsHeadEnd();
939: PetscFunctionReturn(PETSC_SUCCESS);
940: }
942: static PetscErrorCode DMSetUp_Moab(DM dm)
943: {
944: Vec local, global;
945: IS from, to;
946: moab::Range::iterator iter;
947: PetscInt i, j, f, bs, vent, totsize, *lgmap;
948: DM_Moab *dmmoab = (DM_Moab *)dm->data;
949: moab::Range adjs;
951: PetscFunctionBegin;
953: /* Get the local and shared vertices and cache it */
954: PetscCheck(dmmoab->mbiface != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB Interface before calling SetUp.");
955: #if defined(MOAB_HAVE_MPI)
956: PetscCheck(dmmoab->pcomm != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB ParallelComm object before calling SetUp.");
957: #endif
959: /* Get the entities recursively in the current part of the mesh, if user did not set the local vertices explicitly */
960: if (dmmoab->vlocal->empty()) {
961: PetscCallMOAB(dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, 0, *dmmoab->vlocal, false));
963: #if defined(MOAB_HAVE_MPI)
964: /* filter based on parallel status */
965: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned));
967: /* filter all the non-owned and shared entities out of the list */
968: // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
969: adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
970: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost));
971: adjs = moab::subtract(adjs, *dmmoab->vghost);
972: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
973: #else
974: *dmmoab->vowned = *dmmoab->vlocal;
975: #endif
977: /* compute and cache the sizes of local and ghosted entities */
978: dmmoab->nloc = dmmoab->vowned->size();
979: dmmoab->nghost = dmmoab->vghost->size();
981: #if defined(MOAB_HAVE_MPI)
982: PetscCallMPI(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
983: 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));
984: #else
985: dmmoab->n = dmmoab->nloc;
986: #endif
987: }
989: {
990: /* get the information about the local elements in the mesh */
991: dmmoab->eghost->clear();
993: /* first decipher the leading dimension */
994: for (i = 3; i > 0; i--) {
995: dmmoab->elocal->clear();
996: PetscCallMOAB(dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false));
998: /* store the current mesh dimension */
999: if (dmmoab->elocal->size()) {
1000: dmmoab->dim = i;
1001: break;
1002: }
1003: }
1005: PetscCall(DMSetDimension(dm, dmmoab->dim));
1007: #if defined(MOAB_HAVE_MPI)
1008: /* filter the ghosted and owned element list */
1009: *dmmoab->eghost = *dmmoab->elocal;
1010: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT));
1011: *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
1012: #endif
1014: dmmoab->neleloc = dmmoab->elocal->size();
1015: dmmoab->neleghost = dmmoab->eghost->size();
1017: #if defined(MOAB_HAVE_MPI)
1018: PetscCallMPI(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
1019: PetscCall(PetscInfo(NULL, "%d-dim elements: owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost));
1020: #else
1021: dmmoab->nele = dmmoab->neleloc;
1022: #endif
1023: }
1025: bs = dmmoab->bs;
1026: if (!dmmoab->ltog_tag) {
1027: /* Get the global ID tag. The global ID tag is applied to each
1028: vertex. It acts as an global identifier which MOAB uses to
1029: assemble the individual pieces of the mesh */
1030: PetscCallMOAB(dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag));
1031: }
1033: totsize = dmmoab->vlocal->size();
1034: 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);
1035: PetscCall(PetscCalloc1(totsize, &dmmoab->gsindices));
1036: {
1037: /* first get the local indices */
1038: PetscCallMOAB(dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]));
1039: /* next get the ghosted indices */
1040: if (dmmoab->nghost) PetscCallMOAB(dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]));
1042: /* find out the local and global minima of GLOBAL_ID */
1043: dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1044: for (i = 0; i < totsize; ++i) {
1045: if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1046: if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1047: }
1049: PetscCallMPI(MPIU_Allreduce(&dmmoab->lminmax[0], &dmmoab->gminmax[0], 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm));
1050: PetscCallMPI(MPIU_Allreduce(&dmmoab->lminmax[1], &dmmoab->gminmax[1], 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm));
1052: /* set the GID map */
1053: for (i = 0; i < totsize; ++i) dmmoab->gsindices[i] -= dmmoab->gminmax[0]; /* zero based index needed for IS */
1054: dmmoab->lminmax[0] -= dmmoab->gminmax[0];
1055: dmmoab->lminmax[1] -= dmmoab->gminmax[0];
1057: 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]));
1058: }
1059: 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,
1060: dmmoab->numFields);
1062: {
1063: dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1064: dmmoab->seqend = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1065: PetscCall(PetscInfo(NULL, "SEQUENCE: Local [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "]\n", dmmoab->seqstart, dmmoab->seqend));
1067: PetscCall(PetscMalloc2(dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->gidmap, dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->lidmap));
1068: PetscCall(PetscMalloc1(totsize * dmmoab->numFields, &lgmap));
1070: i = j = 0;
1071: /* set the owned vertex data first */
1072: for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1073: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1074: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1075: dmmoab->lidmap[vent] = i;
1076: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1077: }
1078: /* next arrange all the ghosted data information */
1079: for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1080: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1081: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1082: dmmoab->lidmap[vent] = i;
1083: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1084: }
1086: /* We need to create the Global to Local Vector Scatter Contexts
1087: 1) First create a local and global vector
1088: 2) Create a local and global IS
1089: 3) Create VecScatter and LtoGMapping objects
1090: 4) Cleanup the IS and Vec objects
1091: */
1092: PetscCall(DMCreateGlobalVector(dm, &global));
1093: PetscCall(DMCreateLocalVector(dm, &local));
1095: PetscCall(VecGetOwnershipRange(global, &dmmoab->vstart, &dmmoab->vend));
1097: /* global to local must retrieve ghost points */
1098: PetscCall(ISCreateStride(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, dmmoab->vstart, 1, &from));
1099: PetscCall(ISSetBlockSize(from, bs));
1101: PetscCall(ISCreateGeneral(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, &lgmap[0], PETSC_COPY_VALUES, &to));
1102: PetscCall(ISSetBlockSize(to, bs));
1104: if (!dmmoab->ltog_map) {
1105: /* create to the local to global mapping for vectors in order to use VecSetValuesLocal */
1106: PetscCall(ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm, dmmoab->bs, totsize * dmmoab->numFields, lgmap, PETSC_COPY_VALUES, &dmmoab->ltog_map));
1107: }
1109: /* now create the scatter object from local to global vector */
1110: PetscCall(VecScatterCreate(local, from, global, to, &dmmoab->ltog_sendrecv));
1112: /* clean up IS, Vec */
1113: PetscCall(PetscFree(lgmap));
1114: PetscCall(ISDestroy(&from));
1115: PetscCall(ISDestroy(&to));
1116: PetscCall(VecDestroy(&local));
1117: PetscCall(VecDestroy(&global));
1118: }
1120: dmmoab->bndyvtx = new moab::Range();
1121: dmmoab->bndyfaces = new moab::Range();
1122: dmmoab->bndyelems = new moab::Range();
1123: /* skin the boundary and store nodes */
1124: if (!dmmoab->hlevel) {
1125: /* get the skin vertices of boundary faces for the current partition and then filter
1126: the local, boundary faces, vertices and elements alone via PSTATUS flags;
1127: this should not give us any ghosted boundary, but if user needs such a functionality
1128: it would be easy to add it based on the find_skin query below */
1129: moab::Skinner skinner(dmmoab->mbiface);
1131: /* get the entities on the skin - only the faces */
1132: PetscCallMOAB(skinner.find_skin(dmmoab->fileset, *dmmoab->elocal, false, *dmmoab->bndyfaces, NULL, true, true, false));
1134: #if defined(MOAB_HAVE_MPI)
1135: /* filter all the non-owned and shared entities out of the list */
1136: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT));
1137: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT));
1138: #endif
1140: /* get all the nodes via connectivity and the parent elements via adjacency information */
1141: PetscCallMOAB(dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false));
1142: PetscCallMOAB(dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION));
1143: } else {
1144: /* Let us query the hierarchy manager and get the results directly for this level */
1145: for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1146: moab::EntityHandle elemHandle = *iter;
1147: if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1148: dmmoab->bndyelems->insert(elemHandle);
1149: /* For this boundary element, query the vertices and add them to the list */
1150: std::vector<moab::EntityHandle> connect;
1151: PetscCallMOAB(dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect));
1152: for (unsigned iv = 0; iv < connect.size(); ++iv)
1153: if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv])) dmmoab->bndyvtx->insert(connect[iv]);
1154: /* Next, let us query the boundary faces and add them also to the list */
1155: std::vector<moab::EntityHandle> faces;
1156: PetscCallMOAB(dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim - 1, faces));
1157: for (unsigned ifa = 0; ifa < faces.size(); ++ifa)
1158: if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa])) dmmoab->bndyfaces->insert(faces[ifa]);
1159: }
1160: }
1161: #if defined(MOAB_HAVE_MPI)
1162: /* filter all the non-owned and shared entities out of the list */
1163: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT));
1164: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT));
1165: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT));
1166: #endif
1167: }
1168: PetscCall(PetscInfo(NULL, "Found %zu boundary vertices, %zu boundary faces and %zu boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size()));
1170: /* Get the material sets and populate the data for all locally owned elements */
1171: {
1172: PetscCall(PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials));
1173: /* Get the count of entities of particular type from dmmoab->elocal
1174: -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1175: moab::Range msets;
1176: PetscCallMOAB(dmmoab->mbiface->get_entities_by_type_and_tag(dmmoab->fileset, moab::MBENTITYSET, &dmmoab->material_tag, NULL, 1, msets, moab::Interface::UNION));
1177: if (msets.size() == 0) PetscCall(PetscInfo(NULL, "No material sets found in the fileset.\n"));
1179: for (unsigned i = 0; i < msets.size(); ++i) {
1180: moab::Range msetelems;
1181: PetscCallMOAB(dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true));
1182: #if defined(MOAB_HAVE_MPI)
1183: /* filter all the non-owned and shared entities out of the list */
1184: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT));
1185: #endif
1187: int partID;
1188: moab::EntityHandle mset = msets[i];
1189: PetscCallMOAB(dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID));
1191: for (unsigned j = 0; j < msetelems.size(); ++j) dmmoab->materials[dmmoab->elocal->index(msetelems[j])] = partID;
1192: }
1193: }
1194: PetscFunctionReturn(PETSC_SUCCESS);
1195: }
1197: /*@C
1198: DMMoabCreateVertices - Creates and adds several vertices to the primary set represented by the DM.
1200: Collective
1202: Input Parameters:
1203: + dm - The `DM` object
1204: . coords - The connectivity of the element
1205: - nverts - The number of vertices that form the element
1207: Output Parameter:
1208: . overts - The list of vertices that were created (can be `NULL`)
1210: Level: beginner
1212: .seealso: `DMMOAB`, `DMMoabCreateSubmesh()`, `DMMoabCreateElement()`
1213: @*/
1214: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal *coords, PetscInt nverts, moab::Range *overts)
1215: {
1216: DM_Moab *dmmoab;
1217: moab::Range verts;
1219: PetscFunctionBegin;
1221: PetscAssertPointer(coords, 2);
1223: dmmoab = (DM_Moab *)dm->data;
1225: /* Insert new points */
1226: PetscCallMOAB(dmmoab->mbiface->create_vertices(&coords[0], nverts, verts));
1227: PetscCallMOAB(dmmoab->mbiface->add_entities(dmmoab->fileset, verts));
1228: if (overts) *overts = verts;
1229: PetscFunctionReturn(PETSC_SUCCESS);
1230: }
1232: /*@C
1233: DMMoabCreateElement - Adds an element of specified type to the primary set represented by the DM.
1235: Collective
1237: Input Parameters:
1238: + dm - The DM object
1239: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1240: . conn - The connectivity of the element
1241: - nverts - The number of vertices that form the element
1243: Output Parameter:
1244: . oelem - The handle to the element created and added to the `DM` object
1246: Level: beginner
1248: .seealso: `DMMOAB`, `DMMoabCreateSubmesh()`, `DMMoabCreateVertices()`
1249: @*/
1250: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle *conn, PetscInt nverts, moab::EntityHandle *oelem)
1251: {
1252: DM_Moab *dmmoab;
1253: moab::EntityHandle elem;
1255: PetscFunctionBegin;
1257: PetscAssertPointer(conn, 3);
1259: dmmoab = (DM_Moab *)dm->data;
1261: /* Insert new element */
1262: PetscCallMOAB(dmmoab->mbiface->create_element(type, conn, nverts, elem));
1263: PetscCallMOAB(dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1));
1264: if (oelem) *oelem = elem;
1265: PetscFunctionReturn(PETSC_SUCCESS);
1266: }
1268: /*@C
1269: DMMoabCreateSubmesh - Creates a sub-`DM` object with a set that contains all vertices/elements of the parent
1270: in addition to providing support for dynamic mesh modifications. This is useful for AMR calculations to
1271: create a DM object on a refined level.
1273: Collective
1275: Input Parameters:
1276: . dm - The `DM` object
1278: Output Parameter:
1279: . newdm - The sub `DM` object with updated set information
1281: Level: advanced
1283: .seealso: `DMMOAB`, `DMCreate()`, `DMMoabCreateVertices()`, `DMMoabCreateElement()`
1284: @*/
1285: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1286: {
1287: DM_Moab *dmmoab;
1288: DM_Moab *ndmmoab;
1290: PetscFunctionBegin;
1293: dmmoab = (DM_Moab *)dm->data;
1295: /* Create the basic DMMOAB object and keep the default parameters created by DM impls */
1296: PetscCall(DMMoabCreateMoab(((PetscObject)dm)->comm, dmmoab->mbiface, &dmmoab->ltog_tag, NULL, newdm));
1298: /* get all the necessary handles from the private DM object */
1299: ndmmoab = (DM_Moab *)(*newdm)->data;
1301: /* set the sub-mesh's parent DM reference */
1302: ndmmoab->parent = &dm;
1304: /* create a file set to associate all entities in current mesh */
1305: PetscCallMOAB(ndmmoab->mbiface->create_meshset(moab::MESHSET_SET, ndmmoab->fileset));
1307: /* create a meshset and then add old fileset as child */
1308: PetscCallMOAB(ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal));
1309: PetscCallMOAB(ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal));
1311: /* preserve the field association between the parent and sub-mesh objects */
1312: PetscCall(DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames));
1313: PetscFunctionReturn(PETSC_SUCCESS);
1314: }
1316: static PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1317: {
1318: DM_Moab *dmmoab = (DM_Moab *)dm->data;
1319: const char *name;
1320: MPI_Comm comm;
1321: PetscMPIInt size;
1323: PetscFunctionBegin;
1324: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1325: PetscCallMPI(MPI_Comm_size(comm, &size));
1326: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1327: PetscCall(PetscViewerASCIIPushTab(viewer));
1328: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimensions:\n", name, dmmoab->dim));
1329: else PetscCall(PetscViewerASCIIPrintf(viewer, "Mesh in %" PetscInt_FMT " dimensions:\n", dmmoab->dim));
1330: /* print details about the global mesh */
1331: {
1332: PetscCall(PetscViewerASCIIPushTab(viewer));
1333: PetscCall(PetscViewerASCIIPrintf(viewer, "Sizes: cells=%" PetscInt_FMT ", vertices=%" PetscInt_FMT ", blocks=%" PetscInt_FMT "\n", dmmoab->nele, dmmoab->n, dmmoab->bs));
1334: /* print boundary data */
1335: PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary trace:\n"));
1336: {
1337: PetscCall(PetscViewerASCIIPushTab(viewer));
1338: PetscCall(PetscViewerASCIIPrintf(viewer, "cells=%zu, faces=%zu, vertices=%zu\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size()));
1339: PetscCall(PetscViewerASCIIPopTab(viewer));
1340: }
1341: /* print field data */
1342: PetscCall(PetscViewerASCIIPrintf(viewer, "Fields: %" PetscInt_FMT " components\n", dmmoab->numFields));
1343: {
1344: PetscCall(PetscViewerASCIIPushTab(viewer));
1345: for (int i = 0; i < dmmoab->numFields; ++i) PetscCall(PetscViewerASCIIPrintf(viewer, "[%" PetscInt_FMT "] - %s\n", i, dmmoab->fieldNames[i]));
1346: PetscCall(PetscViewerASCIIPopTab(viewer));
1347: }
1348: PetscCall(PetscViewerASCIIPopTab(viewer));
1349: }
1350: PetscCall(PetscViewerASCIIPopTab(viewer));
1351: PetscCall(PetscViewerFlush(viewer));
1352: PetscFunctionReturn(PETSC_SUCCESS);
1353: }
1355: static PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1356: {
1357: PetscFunctionReturn(PETSC_SUCCESS);
1358: }
1360: #if PetscDefined(HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1361: static PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1362: {
1363: PetscFunctionReturn(PETSC_SUCCESS);
1364: }
1365: #endif
1367: static PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1368: {
1369: PetscBool isascii, ishdf5, isvtk;
1371: PetscFunctionBegin;
1374: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1375: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1376: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1377: if (isascii) {
1378: PetscCall(DMMoabView_Ascii(dm, viewer));
1379: } else if (ishdf5) {
1380: #if PetscDefined(HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1381: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ));
1382: PetscCall(DMMoabView_HDF5(dm, viewer));
1383: PetscCall(PetscViewerPopFormat(viewer));
1384: #else
1385: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1386: #endif
1387: } else if (isvtk) PetscCall(DMMoabView_VTK(dm, viewer));
1388: PetscFunctionReturn(PETSC_SUCCESS);
1389: }
1391: static PetscErrorCode DMInitialize_Moab(DM);
1393: static PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1394: {
1395: PetscFunctionBegin;
1396: /* get all the necessary handles from the private DM object */
1397: (*newdm)->data = (DM_Moab *)dm->data;
1398: ((DM_Moab *)dm->data)->refct++;
1400: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMMOAB));
1401: PetscCall(DMInitialize_Moab(*newdm));
1402: PetscFunctionReturn(PETSC_SUCCESS);
1403: }
1405: static PetscErrorCode DMInitialize_Moab(DM dm)
1406: {
1407: PetscFunctionBegin;
1408: dm->ops->view = DMView_Moab;
1409: dm->ops->load = NULL /* DMLoad_Moab */;
1410: dm->ops->setfromoptions = DMSetFromOptions_Moab;
1411: dm->ops->clone = DMClone_Moab;
1412: dm->ops->setup = DMSetUp_Moab;
1413: dm->ops->createlocalsection = NULL;
1414: dm->ops->createsectionpermutation = NULL;
1415: dm->ops->createdefaultconstraints = NULL;
1416: dm->ops->createglobalvector = DMCreateGlobalVector_Moab;
1417: dm->ops->createlocalvector = DMCreateLocalVector_Moab;
1418: dm->ops->getlocaltoglobalmapping = NULL;
1419: dm->ops->createfieldis = NULL;
1420: dm->ops->createcoordinatedm = NULL /* DMCreateCoordinateDM_Moab */;
1421: dm->ops->createcellcoordinatedm = NULL;
1422: dm->ops->getcoloring = NULL;
1423: dm->ops->creatematrix = DMCreateMatrix_Moab;
1424: dm->ops->createinterpolation = DMCreateInterpolation_Moab;
1425: dm->ops->createinjection = NULL /* DMCreateInjection_Moab */;
1426: dm->ops->refine = DMRefine_Moab;
1427: dm->ops->coarsen = DMCoarsen_Moab;
1428: dm->ops->refinehierarchy = DMRefineHierarchy_Moab;
1429: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Moab;
1430: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Moab;
1431: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Moab;
1432: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Moab;
1433: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Moab;
1434: dm->ops->destroy = DMDestroy_Moab;
1435: dm->ops->createsubdm = NULL /* DMCreateSubDM_Moab */;
1436: dm->ops->getdimpoints = NULL /* DMGetDimPoints_Moab */;
1437: dm->ops->locatepoints = NULL /* DMLocatePoints_Moab */;
1438: PetscFunctionReturn(PETSC_SUCCESS);
1439: }
1441: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1442: {
1443: PetscFunctionBegin;
1445: PetscCall(PetscNew((DM_Moab **)&dm->data));
1447: ((DM_Moab *)dm->data)->bs = 1;
1448: ((DM_Moab *)dm->data)->numFields = 1;
1449: ((DM_Moab *)dm->data)->n = 0;
1450: ((DM_Moab *)dm->data)->nloc = 0;
1451: ((DM_Moab *)dm->data)->nghost = 0;
1452: ((DM_Moab *)dm->data)->nele = 0;
1453: ((DM_Moab *)dm->data)->neleloc = 0;
1454: ((DM_Moab *)dm->data)->neleghost = 0;
1455: ((DM_Moab *)dm->data)->ltog_map = NULL;
1456: ((DM_Moab *)dm->data)->ltog_sendrecv = NULL;
1458: ((DM_Moab *)dm->data)->refct = 1;
1459: ((DM_Moab *)dm->data)->parent = NULL;
1460: ((DM_Moab *)dm->data)->vlocal = new moab::Range();
1461: ((DM_Moab *)dm->data)->vowned = new moab::Range();
1462: ((DM_Moab *)dm->data)->vghost = new moab::Range();
1463: ((DM_Moab *)dm->data)->elocal = new moab::Range();
1464: ((DM_Moab *)dm->data)->eghost = new moab::Range();
1466: PetscCall(DMInitialize_Moab(dm));
1467: PetscFunctionReturn(PETSC_SUCCESS);
1468: }