Actual source code: dmmoab.cxx
1: #include <petsc/private/dmmbimpl.h>
3: #include <petscdmmoab.h>
4: #include <MBTagConventions.hpp>
5: #include <moab/NestedRefine.hpp>
6: #include <moab/Skinner.hpp>
8: /*MC
9: DMMOAB = "moab" - A DM object that encapsulates an unstructured mesh described by the MOAB mesh database.
10: Direct access to the MOAB Interface and other mesh manipulation related objects are available
11: through public API. Ability to create global and local representation of Vecs containing all
12: unknowns in the interior and shared boundary via a transparent tag-data wrapper is provided
13: along with utility functions to traverse the mesh and assemble a discrete system via
14: field-based/blocked Vec(Get/Set) methods. Input from and output to different formats are
15: available.
17: Reference: https://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html
19: Level: intermediate
21: .seealso: `DMType`, `DMMoabCreate()`, `DMCreate()`, `DMSetType()`, `DMMoabCreateMoab()`
22: M*/
24: /* External function declarations here */
25: PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
26: PETSC_EXTERN PetscErrorCode DMCreateDefaultConstraints_Moab(DM dm);
27: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
28: PETSC_EXTERN PetscErrorCode DMCreateCoordinateDM_Moab(DM dm, DM *cdm);
29: PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM *dmRefined);
30: PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM *dmCoarsened);
31: PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmRefined[]);
32: PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmCoarsened[]);
33: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm);
34: PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM, Vec *);
35: PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM, Vec *);
36: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
37: PETSC_EXTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM, Vec, InsertMode, Vec);
38: PETSC_EXTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM, Vec, InsertMode, Vec);
39: PETSC_EXTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM, Vec, InsertMode, Vec);
40: PETSC_EXTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM, Vec, InsertMode, Vec);
42: /* Un-implemented routines */
43: /*
44: PETSC_EXTERN PetscErrorCode DMCreatelocalsection_Moab(DM dm);
45: PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dmCoarse, DM dmFine, Mat *mat);
46: PETSC_EXTERN PetscErrorCode DMLoad_Moab(DM dm, PetscViewer viewer);
47: PETSC_EXTERN PetscErrorCode DMGetDimPoints_Moab(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd);
48: PETSC_EXTERN PetscErrorCode DMCreateSubDM_Moab(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
49: PETSC_EXTERN PetscErrorCode DMLocatePoints_Moab(DM dm, Vec v, IS *cellIS);
50: */
52: /*@C
53: DMMoabCreate - Creates a DMMoab object, which encapsulates a moab instance
55: Collective
57: Input Parameter:
58: . comm - The communicator for the DMMoab object
60: Output Parameter:
61: . dmb - The DMMoab object
63: Level: beginner
65: @*/
66: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
67: {
68: PetscFunctionBegin;
70: PetscCall(DMCreate(comm, dmb));
71: PetscCall(DMSetType(*dmb, DMMOAB));
72: PetscFunctionReturn(PETSC_SUCCESS);
73: }
75: /*@C
76: DMMoabCreateMoab - Creates a DMMoab object, optionally from an instance and other data
78: Collective
80: Input Parameters:
81: + comm - The communicator for the DMMoab object
82: . mbiface - (ptr to) the MOAB Instance; if passed in NULL, MOAB instance is created inside PETSc, and destroyed
83: along with the DMMoab
84: . ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
85: - range - If non-NULL, contains range of entities to which DOFs will be assigned
87: Output Parameter:
88: . dmb - The DMMoab object
90: Level: intermediate
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;
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) {
144: PetscCall(DMMoabSetLocalToGlobalTag(dmmb, *ltog_tag));
145: } else {
146: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
147: MBERRNM(merr);
148: if (ltog_tag) *ltog_tag = dmmoab->ltog_tag;
149: }
151: merr = dmmoab->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dmmoab->material_tag);
152: MBERRNM(merr);
154: /* set the local range of entities (vertices) of interest */
155: if (range) PetscCall(DMMoabSetLocalVertices(dmmb, range));
156: *dmb = dmmb;
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: #ifdef MOAB_HAVE_MPI
162: /*@C
163: DMMoabGetParallelComm - Get the ParallelComm used with this DMMoab
165: Collective
167: Input Parameter:
168: . dm - The DMMoab object being set
170: Output Parameter:
171: . pcomm - The ParallelComm for the DMMoab
173: Level: beginner
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: @*/
198: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
199: {
200: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
202: PetscFunctionBegin;
205: #ifdef MOAB_HAVE_MPI
206: dmmoab->pcomm = NULL;
207: #endif
208: dmmoab->mbiface = mbiface;
209: dmmoab->icreatedinstance = PETSC_FALSE;
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: /*@C
214: DMMoabGetInterface - Get the MOAB instance used with this DMMoab
216: Collective
218: Input Parameter:
219: . dm - The DMMoab object being set
221: Output Parameter:
222: . mbiface - The MOAB instance set on this DMMoab
224: Level: beginner
226: @*/
227: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
228: {
229: static PetscBool cite = PETSC_FALSE;
231: PetscFunctionBegin;
233: 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, "
234: "K. and Stimpson, C. and Ernst, C.},\n year = {2004}, note = {Report}\n}\n",
235: &cite));
236: *mbiface = ((DM_Moab *)dm->data)->mbiface;
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /*@C
241: DMMoabSetLocalVertices - Set the entities having DOFs on this DMMoab
243: Collective
245: Input Parameters:
246: + dm - The DMMoab object being set
247: - range - The entities treated by this DMMoab
249: Level: beginner
251: @*/
252: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
253: {
254: moab::Range tmpvtxs;
255: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
257: PetscFunctionBegin;
259: dmmoab->vlocal->clear();
260: dmmoab->vowned->clear();
262: dmmoab->vlocal->insert(range->begin(), range->end());
264: #ifdef MOAB_HAVE_MPI
265: moab::ErrorCode merr;
266: /* filter based on parallel status */
267: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
268: MBERRNM(merr);
270: /* filter all the non-owned and shared entities out of the list */
271: tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
272: merr = dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
273: MBERRNM(merr);
274: tmpvtxs = moab::subtract(tmpvtxs, *dmmoab->vghost);
275: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
276: #else
277: *dmmoab->vowned = *dmmoab->vlocal;
278: #endif
280: /* compute and cache the sizes of local and ghosted entities */
281: dmmoab->nloc = dmmoab->vowned->size();
282: dmmoab->nghost = dmmoab->vghost->size();
283: #ifdef MOAB_HAVE_MPI
284: PetscCall(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
285: #else
286: dmmoab->n = dmmoab->nloc;
287: #endif
288: PetscFunctionReturn(PETSC_SUCCESS);
289: }
291: /*@C
292: DMMoabGetAllVertices - Get the entities having DOFs on this DMMoab
294: Collective
296: Input Parameter:
297: . dm - The DMMoab object being set
299: Output Parameter:
300: . owned - The local vertex entities in this DMMoab = (owned+ghosted)
302: Level: beginner
304: @*/
305: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
306: {
307: PetscFunctionBegin;
309: if (local) *local = *((DM_Moab *)dm->data)->vlocal;
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: /*@C
314: DMMoabGetLocalVertices - Get the entities having DOFs on this DMMoab
316: Collective
318: Input Parameter:
319: . dm - The DMMoab object being set
321: Output Parameters:
322: + owned - The owned vertex entities in this DMMoab
323: - ghost - The ghosted entities (non-owned) stored locally in this partition
325: Level: beginner
327: @*/
328: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
329: {
330: PetscFunctionBegin;
332: if (owned) *owned = ((DM_Moab *)dm->data)->vowned;
333: if (ghost) *ghost = ((DM_Moab *)dm->data)->vghost;
334: PetscFunctionReturn(PETSC_SUCCESS);
335: }
337: /*@C
338: DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned
340: Collective
342: Input Parameter:
343: . dm - The DMMoab object being set
345: Output Parameter:
346: . range - The entities owned locally
348: Level: beginner
350: @*/
351: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
352: {
353: PetscFunctionBegin;
355: if (range) *range = ((DM_Moab *)dm->data)->elocal;
356: PetscFunctionReturn(PETSC_SUCCESS);
357: }
359: /*@C
360: DMMoabSetLocalElements - Set the entities having DOFs on this DMMoab
362: Collective
364: Input Parameters:
365: + dm - The DMMoab object being set
366: - range - The entities treated by this DMMoab
368: Level: beginner
370: @*/
371: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
372: {
373: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
375: PetscFunctionBegin;
377: dmmoab->elocal->clear();
378: dmmoab->eghost->clear();
379: dmmoab->elocal->insert(range->begin(), range->end());
380: #ifdef MOAB_HAVE_MPI
381: moab::ErrorCode merr;
382: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
383: MBERRNM(merr);
384: *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
385: #endif
386: dmmoab->neleloc = dmmoab->elocal->size();
387: dmmoab->neleghost = dmmoab->eghost->size();
388: #ifdef MOAB_HAVE_MPI
389: PetscCall(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
390: PetscCall(PetscInfo(dm, "Created %" PetscInt_FMT " local and %" PetscInt_FMT " global elements.\n", dmmoab->neleloc, dmmoab->nele));
391: #else
392: dmmoab->nele = dmmoab->neleloc;
393: #endif
394: PetscFunctionReturn(PETSC_SUCCESS);
395: }
397: /*@C
398: DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
400: Collective
402: Input Parameters:
403: + dm - The DMMoab object being set
404: - ltogtag - The MOAB tag used for local to global ids
406: Level: beginner
408: @*/
409: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
410: {
411: PetscFunctionBegin;
413: ((DM_Moab *)dm->data)->ltog_tag = ltogtag;
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: /*@C
418: DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
420: Collective
422: Input Parameter:
423: . dm - The DMMoab object being set
425: Output Parameter:
426: . ltogtag - The MOAB tag used for local to global ids
428: Level: beginner
430: @*/
431: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
432: {
433: PetscFunctionBegin;
435: *ltog_tag = ((DM_Moab *)dm->data)->ltog_tag;
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }
439: /*@C
440: DMMoabSetBlockSize - Set the block size used with this DMMoab
442: Collective
444: Input Parameters:
445: + dm - The DMMoab object being set
446: - bs - The block size used with this DMMoab
448: Level: beginner
450: @*/
451: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
452: {
453: PetscFunctionBegin;
455: ((DM_Moab *)dm->data)->bs = bs;
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
459: /*@C
460: DMMoabGetBlockSize - Get the block size used with this DMMoab
462: Collective
464: Input Parameter:
465: . dm - The DMMoab object being set
467: Output Parameter:
468: . bs - The block size used with this DMMoab
470: Level: beginner
472: @*/
473: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
474: {
475: PetscFunctionBegin;
477: *bs = ((DM_Moab *)dm->data)->bs;
478: PetscFunctionReturn(PETSC_SUCCESS);
479: }
481: /*@C
482: DMMoabGetSize - Get the global vertex size used with this DMMoab
484: Collective on dm
486: Input Parameter:
487: . dm - The DMMoab object being set
489: Output Parameters:
490: + neg - The number of global elements in the DMMoab instance
491: - nvg - The number of global vertices in the DMMoab instance
493: Level: beginner
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 on dm
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: @*/
522: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
523: {
524: PetscFunctionBegin;
526: if (nel) *nel = ((DM_Moab *)dm->data)->neleloc;
527: if (neg) *neg = ((DM_Moab *)dm->data)->neleghost;
528: if (nvl) *nvl = ((DM_Moab *)dm->data)->nloc;
529: if (nvg) *nvg = ((DM_Moab *)dm->data)->nghost;
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: /*@C
534: DMMoabGetOffset - Get the local offset for the global vector
536: Collective
538: Input Parameter:
539: . dm - The DMMoab object being set
541: Output Parameter:
542: . offset - The local offset for the global vector
544: Level: beginner
546: @*/
547: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
548: {
549: PetscFunctionBegin;
551: *offset = ((DM_Moab *)dm->data)->vstart;
552: PetscFunctionReturn(PETSC_SUCCESS);
553: }
555: /*@C
556: DMMoabGetDimension - Get the dimension of the DM Mesh
558: Collective
560: Input Parameter:
561: . dm - The DMMoab object
563: Output Parameter:
564: . dim - The dimension of DM
566: Level: beginner
568: @*/
569: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
570: {
571: PetscFunctionBegin;
573: *dim = ((DM_Moab *)dm->data)->dim;
574: PetscFunctionReturn(PETSC_SUCCESS);
575: }
577: /*@C
578: DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
579: generated through uniform refinement.
581: Collective on dm
583: Input Parameter:
584: . dm - The DMMoab object being set
586: Output Parameter:
587: . nvg - The current mesh hierarchy level
589: Level: beginner
591: @*/
592: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
593: {
594: PetscFunctionBegin;
596: if (nlevel) *nlevel = ((DM_Moab *)dm->data)->hlevel;
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: /*@C
601: DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh
603: Collective
605: Input Parameters:
606: + dm - The DMMoab object
607: - ehandle - The element entity handle
609: Output Parameter:
610: . mat - The material ID for the current entity
612: Level: beginner
614: @*/
615: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
616: {
617: DM_Moab *dmmoab;
619: PetscFunctionBegin;
621: if (*mat) {
622: dmmoab = (DM_Moab *)(dm)->data;
623: *mat = dmmoab->materials[dmmoab->elocal->index(ehandle)];
624: }
625: PetscFunctionReturn(PETSC_SUCCESS);
626: }
628: /*@C
629: DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities
631: Collective
633: Input Parameters:
634: + dm - The DMMoab object
635: . nconn - Number of entities whose coordinates are needed
636: - conn - The vertex entity handles
638: Output Parameter:
639: . vpos - The coordinates of the requested vertex entities
641: Level: beginner
643: .seealso: `DMMoabGetVertexConnectivity()`
644: @*/
645: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
646: {
647: DM_Moab *dmmoab;
648: moab::ErrorCode merr;
650: PetscFunctionBegin;
654: dmmoab = (DM_Moab *)(dm)->data;
656: /* Get connectivity information in MOAB canonical ordering */
657: if (dmmoab->hlevel) {
658: merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle *>(conn), nconn, dmmoab->hlevel, vpos);
659: MBERRNM(merr);
660: } else {
661: merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);
662: MBERRNM(merr);
663: }
664: PetscFunctionReturn(PETSC_SUCCESS);
665: }
667: /*@C
668: DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity
670: Collective
672: Input Parameters:
673: + dm - The DMMoab object
674: - vhandle - Vertex entity handle
676: Output Parameters:
677: + nconn - Number of entities whose coordinates are needed
678: - conn - The vertex entity handles
680: Level: beginner
682: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabRestoreVertexConnectivity()`
683: @*/
684: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt *nconn, moab::EntityHandle **conn)
685: {
686: DM_Moab *dmmoab;
687: std::vector<moab::EntityHandle> adj_entities, connect;
688: moab::ErrorCode merr;
690: PetscFunctionBegin;
693: dmmoab = (DM_Moab *)(dm)->data;
695: /* Get connectivity information in MOAB canonical ordering */
696: merr = dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION);
697: MBERRNM(merr);
698: merr = dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect);
699: MBERRNM(merr);
701: if (conn) {
702: PetscCall(PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn));
703: PetscCall(PetscArraycpy(*conn, &connect[0], connect.size()));
704: }
705: if (nconn) *nconn = connect.size();
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: /*@C
710: DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity
712: Collective
714: Input Parameters:
715: + dm - The DMMoab object
716: . vhandle - Vertex entity handle
717: . nconn - Number of entities whose coordinates are needed
718: - conn - The vertex entity handles
720: Level: beginner
722: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`
723: @*/
724: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, moab::EntityHandle **conn)
725: {
726: PetscFunctionBegin;
730: if (conn) PetscCall(PetscFree(*conn));
731: if (nconn) *nconn = 0;
732: PetscFunctionReturn(PETSC_SUCCESS);
733: }
735: /*@C
736: DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity
738: Collective
740: Input Parameters:
741: + dm - The DMMoab object
742: - ehandle - Vertex entity handle
744: Output Parameters:
745: + nconn - Number of entities whose coordinates are needed
746: - conn - The vertex entity handles
748: Level: beginner
750: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`, `DMMoabRestoreVertexConnectivity()`
751: @*/
752: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, const moab::EntityHandle **conn)
753: {
754: DM_Moab *dmmoab;
755: const moab::EntityHandle *connect;
756: std::vector<moab::EntityHandle> vconn;
757: moab::ErrorCode merr;
758: PetscInt nnodes;
760: PetscFunctionBegin;
763: dmmoab = (DM_Moab *)(dm)->data;
765: /* Get connectivity information in MOAB canonical ordering */
766: merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes);
767: MBERRNM(merr);
768: if (conn) *conn = connect;
769: if (nconn) *nconn = nnodes;
770: PetscFunctionReturn(PETSC_SUCCESS);
771: }
773: /*@C
774: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
776: Collective
778: Input Parameters:
779: + dm - The DMMoab object
780: - ent - Entity handle
782: Output Parameter:
783: . ent_on_boundary - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
785: Level: beginner
787: .seealso: `DMMoabCheckBoundaryVertices()`
788: @*/
789: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool *ent_on_boundary)
790: {
791: moab::EntityType etype;
792: DM_Moab *dmmoab;
793: PetscInt edim;
795: PetscFunctionBegin;
798: dmmoab = (DM_Moab *)(dm)->data;
800: /* get the entity type and handle accordingly */
801: etype = dmmoab->mbiface->type_from_handle(ent);
802: PetscCheck(etype < moab::MBPOLYHEDRON, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Entity type on the boundary skin is invalid. EntityType = %" PetscInt_FMT, etype);
804: /* get the entity dimension */
805: edim = dmmoab->mbiface->dimension_from_handle(ent);
807: *ent_on_boundary = PETSC_FALSE;
808: if (etype == moab::MBVERTEX && edim == 0) {
809: *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
810: } else {
811: if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
812: if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
813: } else { /* next check the lower-dimensional faces */
814: if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
815: }
816: }
817: PetscFunctionReturn(PETSC_SUCCESS);
818: }
820: /*@C
821: DMMoabCheckBoundaryVertices - Check whether a given entity is on the boundary (vertex, edge, face, element)
823: Input Parameters:
824: + dm - The DMMoab object
825: . nconn - Number of handles
826: - cnt - Array of entity handles
828: Output Parameter:
829: . isbdvtx - Array of boundary markers - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
831: Level: beginner
833: .seealso: `DMMoabIsEntityOnBoundary()`
834: @*/
835: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool *isbdvtx)
836: {
837: DM_Moab *dmmoab;
838: PetscInt i;
840: PetscFunctionBegin;
844: dmmoab = (DM_Moab *)(dm)->data;
846: for (i = 0; i < nconn; ++i) isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
847: PetscFunctionReturn(PETSC_SUCCESS);
848: }
850: /*@C
851: DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary
853: Input Parameter:
854: . dm - The DMMoab object
856: Output Parameters:
857: + bdvtx - Boundary vertices
858: . bdelems - Boundary elements
859: - bdfaces - Boundary faces
861: Level: beginner
863: .seealso: `DMMoabCheckBoundaryVertices()`, `DMMoabIsEntityOnBoundary()`
864: @*/
865: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range **bdelems, const moab::Range **bdfaces)
866: {
867: DM_Moab *dmmoab;
869: PetscFunctionBegin;
871: dmmoab = (DM_Moab *)(dm)->data;
873: if (bdvtx) *bdvtx = dmmoab->bndyvtx;
874: if (bdfaces) *bdfaces = dmmoab->bndyfaces;
875: if (bdelems) *bdfaces = dmmoab->bndyelems;
876: PetscFunctionReturn(PETSC_SUCCESS);
877: }
879: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
880: {
881: PetscInt i;
882: moab::ErrorCode merr;
883: DM_Moab *dmmoab = (DM_Moab *)dm->data;
885: PetscFunctionBegin;
888: dmmoab->refct--;
889: if (!dmmoab->refct) {
890: delete dmmoab->vlocal;
891: delete dmmoab->vowned;
892: delete dmmoab->vghost;
893: delete dmmoab->elocal;
894: delete dmmoab->eghost;
895: delete dmmoab->bndyvtx;
896: delete dmmoab->bndyfaces;
897: delete dmmoab->bndyelems;
899: PetscCall(PetscFree(dmmoab->gsindices));
900: PetscCall(PetscFree2(dmmoab->gidmap, dmmoab->lidmap));
901: PetscCall(PetscFree(dmmoab->dfill));
902: PetscCall(PetscFree(dmmoab->ofill));
903: PetscCall(PetscFree(dmmoab->materials));
904: if (dmmoab->fieldNames) {
905: for (i = 0; i < dmmoab->numFields; i++) PetscCall(PetscFree(dmmoab->fieldNames[i]));
906: PetscCall(PetscFree(dmmoab->fieldNames));
907: }
909: if (dmmoab->nhlevels) {
910: PetscCall(PetscFree(dmmoab->hsets));
911: dmmoab->nhlevels = 0;
912: if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
913: dmmoab->hierarchy = NULL;
914: }
916: if (dmmoab->icreatedinstance) {
917: delete dmmoab->pcomm;
918: merr = dmmoab->mbiface->delete_mesh();
919: MBERRNM(merr);
920: delete dmmoab->mbiface;
921: }
922: dmmoab->mbiface = NULL;
923: #ifdef MOAB_HAVE_MPI
924: dmmoab->pcomm = NULL;
925: #endif
926: PetscCall(VecScatterDestroy(&dmmoab->ltog_sendrecv));
927: PetscCall(ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map));
928: PetscCall(PetscFree(dm->data));
929: }
930: PetscFunctionReturn(PETSC_SUCCESS);
931: }
933: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(DM dm, PetscOptionItems *PetscOptionsObject)
934: {
935: DM_Moab *dmmoab = (DM_Moab *)dm->data;
937: PetscFunctionBegin;
938: PetscOptionsHeadBegin(PetscOptionsObject, "DMMoab Options");
939: PetscCall(PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL, 0));
940: 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));
941: /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
942: 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));
943: 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));
944: PetscCall(PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum *)&dmmoab->read_mode, NULL));
945: PetscCall(PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum *)&dmmoab->write_mode, NULL));
946: PetscOptionsHeadEnd();
947: PetscFunctionReturn(PETSC_SUCCESS);
948: }
950: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
951: {
952: moab::ErrorCode merr;
953: Vec local, global;
954: IS from, to;
955: moab::Range::iterator iter;
956: PetscInt i, j, f, bs, vent, totsize, *lgmap;
957: DM_Moab *dmmoab = (DM_Moab *)dm->data;
958: moab::Range adjs;
960: PetscFunctionBegin;
962: /* Get the local and shared vertices and cache it */
963: PetscCheck(dmmoab->mbiface != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB Interface before calling SetUp.");
964: #ifdef MOAB_HAVE_MPI
965: PetscCheck(dmmoab->pcomm != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB ParallelComm object before calling SetUp.");
966: #endif
968: /* Get the entities recursively in the current part of the mesh, if user did not set the local vertices explicitly */
969: if (dmmoab->vlocal->empty()) {
970: //merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,*dmmoab->vlocal,true);MBERRNM(merr);
971: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, 0, *dmmoab->vlocal, false);
972: MBERRNM(merr);
974: #ifdef MOAB_HAVE_MPI
975: /* filter based on parallel status */
976: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
977: MBERRNM(merr);
979: /* filter all the non-owned and shared entities out of the list */
980: // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
981: adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
982: merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
983: MBERRNM(merr);
984: adjs = moab::subtract(adjs, *dmmoab->vghost);
985: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
986: #else
987: *dmmoab->vowned = *dmmoab->vlocal;
988: #endif
990: /* compute and cache the sizes of local and ghosted entities */
991: dmmoab->nloc = dmmoab->vowned->size();
992: dmmoab->nghost = dmmoab->vghost->size();
994: #ifdef MOAB_HAVE_MPI
995: PetscCall(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
996: 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));
997: #else
998: dmmoab->n = dmmoab->nloc;
999: #endif
1000: }
1002: {
1003: /* get the information about the local elements in the mesh */
1004: dmmoab->eghost->clear();
1006: /* first decipher the leading dimension */
1007: for (i = 3; i > 0; i--) {
1008: dmmoab->elocal->clear();
1009: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false);
1010: MBERRNM(merr);
1012: /* store the current mesh dimension */
1013: if (dmmoab->elocal->size()) {
1014: dmmoab->dim = i;
1015: break;
1016: }
1017: }
1019: PetscCall(DMSetDimension(dm, dmmoab->dim));
1021: #ifdef MOAB_HAVE_MPI
1022: /* filter the ghosted and owned element list */
1023: *dmmoab->eghost = *dmmoab->elocal;
1024: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1025: MBERRNM(merr);
1026: *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
1027: #endif
1029: dmmoab->neleloc = dmmoab->elocal->size();
1030: dmmoab->neleghost = dmmoab->eghost->size();
1032: #ifdef MOAB_HAVE_MPI
1033: PetscCall(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
1034: PetscCall(PetscInfo(NULL, "%d-dim elements: owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost));
1035: #else
1036: dmmoab->nele = dmmoab->neleloc;
1037: #endif
1038: }
1040: bs = dmmoab->bs;
1041: if (!dmmoab->ltog_tag) {
1042: /* Get the global ID tag. The global ID tag is applied to each
1043: vertex. It acts as an global identifier which MOAB uses to
1044: assemble the individual pieces of the mesh */
1045: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
1046: MBERRNM(merr);
1047: }
1049: totsize = dmmoab->vlocal->size();
1050: 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);
1051: PetscCall(PetscCalloc1(totsize, &dmmoab->gsindices));
1052: {
1053: /* first get the local indices */
1054: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]);
1055: MBERRNM(merr);
1056: if (dmmoab->nghost) { /* next get the ghosted indices */
1057: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]);
1058: MBERRNM(merr);
1059: }
1061: /* find out the local and global minima of GLOBAL_ID */
1062: dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1063: for (i = 0; i < totsize; ++i) {
1064: if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1065: if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1066: }
1068: PetscCall(MPIU_Allreduce(&dmmoab->lminmax[0], &dmmoab->gminmax[0], 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm));
1069: PetscCall(MPIU_Allreduce(&dmmoab->lminmax[1], &dmmoab->gminmax[1], 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm));
1071: /* set the GID map */
1072: for (i = 0; i < totsize; ++i) { dmmoab->gsindices[i] -= dmmoab->gminmax[0]; /* zero based index needed for IS */ }
1073: dmmoab->lminmax[0] -= dmmoab->gminmax[0];
1074: dmmoab->lminmax[1] -= dmmoab->gminmax[0];
1076: 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]));
1077: }
1078: 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,
1079: dmmoab->numFields);
1081: {
1082: dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1083: dmmoab->seqend = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1084: PetscCall(PetscInfo(NULL, "SEQUENCE: Local [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "]\n", dmmoab->seqstart, dmmoab->seqend));
1086: PetscCall(PetscMalloc2(dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->gidmap, dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->lidmap));
1087: PetscCall(PetscMalloc1(totsize * dmmoab->numFields, &lgmap));
1089: i = j = 0;
1090: /* set the owned vertex data first */
1091: for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1092: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1093: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1094: dmmoab->lidmap[vent] = i;
1095: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1096: }
1097: /* next arrange all the ghosted data information */
1098: for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1099: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1100: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1101: dmmoab->lidmap[vent] = i;
1102: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1103: }
1105: /* We need to create the Global to Local Vector Scatter Contexts
1106: 1) First create a local and global vector
1107: 2) Create a local and global IS
1108: 3) Create VecScatter and LtoGMapping objects
1109: 4) Cleanup the IS and Vec objects
1110: */
1111: PetscCall(DMCreateGlobalVector(dm, &global));
1112: PetscCall(DMCreateLocalVector(dm, &local));
1114: PetscCall(VecGetOwnershipRange(global, &dmmoab->vstart, &dmmoab->vend));
1116: /* global to local must retrieve ghost points */
1117: PetscCall(ISCreateStride(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, dmmoab->vstart, 1, &from));
1118: PetscCall(ISSetBlockSize(from, bs));
1120: PetscCall(ISCreateGeneral(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, &lgmap[0], PETSC_COPY_VALUES, &to));
1121: PetscCall(ISSetBlockSize(to, bs));
1123: if (!dmmoab->ltog_map) {
1124: /* create to the local to global mapping for vectors in order to use VecSetValuesLocal */
1125: PetscCall(ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm, dmmoab->bs, totsize * dmmoab->numFields, lgmap, PETSC_COPY_VALUES, &dmmoab->ltog_map));
1126: }
1128: /* now create the scatter object from local to global vector */
1129: PetscCall(VecScatterCreate(local, from, global, to, &dmmoab->ltog_sendrecv));
1131: /* clean up IS, Vec */
1132: PetscCall(PetscFree(lgmap));
1133: PetscCall(ISDestroy(&from));
1134: PetscCall(ISDestroy(&to));
1135: PetscCall(VecDestroy(&local));
1136: PetscCall(VecDestroy(&global));
1137: }
1139: dmmoab->bndyvtx = new moab::Range();
1140: dmmoab->bndyfaces = new moab::Range();
1141: dmmoab->bndyelems = new moab::Range();
1142: /* skin the boundary and store nodes */
1143: if (!dmmoab->hlevel) {
1144: /* get the skin vertices of boundary faces for the current partition and then filter
1145: the local, boundary faces, vertices and elements alone via PSTATUS flags;
1146: this should not give us any ghosted boundary, but if user needs such a functionality
1147: it would be easy to add it based on the find_skin query below */
1148: moab::Skinner skinner(dmmoab->mbiface);
1150: /* get the entities on the skin - only the faces */
1151: merr = skinner.find_skin(dmmoab->fileset, *dmmoab->elocal, false, *dmmoab->bndyfaces, NULL, true, true, false);
1152: MBERRNM(merr); // 'false' param indicates we want faces back, not vertices
1154: #ifdef MOAB_HAVE_MPI
1155: /* filter all the non-owned and shared entities out of the list */
1156: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1157: MBERRNM(merr);
1158: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT);
1159: MBERRNM(merr);
1160: #endif
1162: /* get all the nodes via connectivity and the parent elements via adjacency information */
1163: merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false);
1164: MBERRNM(merr);
1165: merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION);
1166: MBERRNM(merr);
1167: } else {
1168: /* Let us query the hierarchy manager and get the results directly for this level */
1169: for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1170: moab::EntityHandle elemHandle = *iter;
1171: if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1172: dmmoab->bndyelems->insert(elemHandle);
1173: /* For this boundary element, query the vertices and add them to the list */
1174: std::vector<moab::EntityHandle> connect;
1175: merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect);
1176: MBERRNM(merr);
1177: for (unsigned iv = 0; iv < connect.size(); ++iv)
1178: if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv])) dmmoab->bndyvtx->insert(connect[iv]);
1179: /* Next, let us query the boundary faces and add them also to the list */
1180: std::vector<moab::EntityHandle> faces;
1181: merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim - 1, faces);
1182: MBERRNM(merr);
1183: for (unsigned ifa = 0; ifa < faces.size(); ++ifa)
1184: if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa])) dmmoab->bndyfaces->insert(faces[ifa]);
1185: }
1186: }
1187: #ifdef MOAB_HAVE_MPI
1188: /* filter all the non-owned and shared entities out of the list */
1189: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1190: MBERRNM(merr);
1191: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1192: MBERRNM(merr);
1193: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1194: MBERRNM(merr);
1195: #endif
1196: }
1197: PetscCall(PetscInfo(NULL, "Found %zu boundary vertices, %zu boundary faces and %zu boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size()));
1199: /* Get the material sets and populate the data for all locally owned elements */
1200: {
1201: PetscCall(PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials));
1202: /* Get the count of entities of particular type from dmmoab->elocal
1203: -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1204: moab::Range msets;
1205: merr = dmmoab->mbiface->get_entities_by_type_and_tag(dmmoab->fileset, moab::MBENTITYSET, &dmmoab->material_tag, NULL, 1, msets, moab::Interface::UNION);
1206: MBERRNM(merr);
1207: if (msets.size() == 0) PetscCall(PetscInfo(NULL, "No material sets found in the fileset."));
1209: for (unsigned i = 0; i < msets.size(); ++i) {
1210: moab::Range msetelems;
1211: merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);
1212: MBERRNM(merr);
1213: #ifdef MOAB_HAVE_MPI
1214: /* filter all the non-owned and shared entities out of the list */
1215: merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1216: MBERRNM(merr);
1217: #endif
1219: int partID;
1220: moab::EntityHandle mset = msets[i];
1221: merr = dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID);
1222: MBERRNM(merr);
1224: for (unsigned j = 0; j < msetelems.size(); ++j) dmmoab->materials[dmmoab->elocal->index(msetelems[j])] = partID;
1225: }
1226: }
1228: PetscFunctionReturn(PETSC_SUCCESS);
1229: }
1231: /*@C
1232: DMMoabCreateVertices - Creates and adds several vertices to the primary set represented by the DM.
1234: Collective
1236: Input Parameters:
1237: + dm - The DM object
1238: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1239: . conn - The connectivity of the element
1240: - nverts - The number of vertices that form the element
1242: Output Parameter:
1243: . overts - The list of vertices that were created (can be NULL)
1245: Level: beginner
1247: .seealso: `DMMoabCreateSubmesh()`, `DMMoabCreateElement()`
1248: @*/
1249: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal *coords, PetscInt nverts, moab::Range *overts)
1250: {
1251: moab::ErrorCode merr;
1252: DM_Moab *dmmoab;
1253: moab::Range verts;
1255: PetscFunctionBegin;
1259: dmmoab = (DM_Moab *)dm->data;
1261: /* Insert new points */
1262: merr = dmmoab->mbiface->create_vertices(&coords[0], nverts, verts);
1263: MBERRNM(merr);
1264: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, verts);
1265: MBERRNM(merr);
1267: if (overts) *overts = verts;
1268: PetscFunctionReturn(PETSC_SUCCESS);
1269: }
1271: /*@C
1272: DMMoabCreateElement - Adds an element of specified type to the primary set represented by the DM.
1274: Collective
1276: Input Parameters:
1277: + dm - The DM object
1278: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1279: . conn - The connectivity of the element
1280: - nverts - The number of vertices that form the element
1282: Output Parameter:
1283: . oelem - The handle to the element created and added to the DM object
1285: Level: beginner
1287: .seealso: `DMMoabCreateSubmesh()`, `DMMoabCreateVertices()`
1288: @*/
1289: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle *conn, PetscInt nverts, moab::EntityHandle *oelem)
1290: {
1291: moab::ErrorCode merr;
1292: DM_Moab *dmmoab;
1293: moab::EntityHandle elem;
1295: PetscFunctionBegin;
1299: dmmoab = (DM_Moab *)dm->data;
1301: /* Insert new element */
1302: merr = dmmoab->mbiface->create_element(type, conn, nverts, elem);
1303: MBERRNM(merr);
1304: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1);
1305: MBERRNM(merr);
1307: if (oelem) *oelem = elem;
1308: PetscFunctionReturn(PETSC_SUCCESS);
1309: }
1311: /*@C
1312: DMMoabCreateSubmesh - Creates a sub-DM object with a set that contains all vertices/elements of the parent
1313: in addition to providing support for dynamic mesh modifications. This is useful for AMR calculations to
1314: create a DM object on a refined level.
1316: Collective
1318: Input Parameters:
1319: . dm - The DM object
1321: Output Parameter:
1322: . newdm - The sub DM object with updated set information
1324: Level: advanced
1326: .seealso: `DMCreate()`, `DMMoabCreateVertices()`, `DMMoabCreateElement()`
1327: @*/
1328: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1329: {
1330: DM_Moab *dmmoab;
1331: DM_Moab *ndmmoab;
1332: moab::ErrorCode merr;
1334: PetscFunctionBegin;
1337: dmmoab = (DM_Moab *)dm->data;
1339: /* Create the basic DMMoab object and keep the default parameters created by DM impls */
1340: PetscCall(DMMoabCreateMoab(((PetscObject)dm)->comm, dmmoab->mbiface, &dmmoab->ltog_tag, NULL, newdm));
1342: /* get all the necessary handles from the private DM object */
1343: ndmmoab = (DM_Moab *)(*newdm)->data;
1345: /* set the sub-mesh's parent DM reference */
1346: ndmmoab->parent = &dm;
1348: /* create a file set to associate all entities in current mesh */
1349: merr = ndmmoab->mbiface->create_meshset(moab::MESHSET_SET, ndmmoab->fileset);
1350: MBERR("Creating file set failed", merr);
1352: /* create a meshset and then add old fileset as child */
1353: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal);
1354: MBERR("Adding child vertices to parent failed", merr);
1355: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal);
1356: MBERR("Adding child elements to parent failed", merr);
1358: /* preserve the field association between the parent and sub-mesh objects */
1359: PetscCall(DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames));
1360: PetscFunctionReturn(PETSC_SUCCESS);
1361: }
1363: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1364: {
1365: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
1366: const char *name;
1367: MPI_Comm comm;
1368: PetscMPIInt size;
1370: PetscFunctionBegin;
1371: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1372: PetscCallMPI(MPI_Comm_size(comm, &size));
1373: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1374: PetscCall(PetscViewerASCIIPushTab(viewer));
1375: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimensions:\n", name, dmmoab->dim));
1376: else PetscCall(PetscViewerASCIIPrintf(viewer, "Mesh in %" PetscInt_FMT " dimensions:\n", dmmoab->dim));
1377: /* print details about the global mesh */
1378: {
1379: PetscCall(PetscViewerASCIIPushTab(viewer));
1380: PetscCall(PetscViewerASCIIPrintf(viewer, "Sizes: cells=%" PetscInt_FMT ", vertices=%" PetscInt_FMT ", blocks=%" PetscInt_FMT "\n", dmmoab->nele, dmmoab->n, dmmoab->bs));
1381: /* print boundary data */
1382: PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary trace:\n"));
1383: {
1384: PetscCall(PetscViewerASCIIPushTab(viewer));
1385: PetscCall(PetscViewerASCIIPrintf(viewer, "cells=%zu, faces=%zu, vertices=%zu\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size()));
1386: PetscCall(PetscViewerASCIIPopTab(viewer));
1387: }
1388: /* print field data */
1389: PetscCall(PetscViewerASCIIPrintf(viewer, "Fields: %" PetscInt_FMT " components\n", dmmoab->numFields));
1390: {
1391: PetscCall(PetscViewerASCIIPushTab(viewer));
1392: for (int i = 0; i < dmmoab->numFields; ++i) PetscCall(PetscViewerASCIIPrintf(viewer, "[%" PetscInt_FMT "] - %s\n", i, dmmoab->fieldNames[i]));
1393: PetscCall(PetscViewerASCIIPopTab(viewer));
1394: }
1395: PetscCall(PetscViewerASCIIPopTab(viewer));
1396: }
1397: PetscCall(PetscViewerASCIIPopTab(viewer));
1398: PetscCall(PetscViewerFlush(viewer));
1399: PetscFunctionReturn(PETSC_SUCCESS);
1400: }
1402: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1403: {
1404: PetscFunctionReturn(PETSC_SUCCESS);
1405: }
1407: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1408: {
1409: PetscFunctionReturn(PETSC_SUCCESS);
1410: }
1412: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1413: {
1414: PetscBool iascii, ishdf5, isvtk;
1416: PetscFunctionBegin;
1419: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1420: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1421: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1422: if (iascii) {
1423: PetscCall(DMMoabView_Ascii(dm, viewer));
1424: } else if (ishdf5) {
1425: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1426: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ));
1427: PetscCall(DMMoabView_HDF5(dm, viewer));
1428: PetscCall(PetscViewerPopFormat(viewer));
1429: #else
1430: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1431: #endif
1432: } else if (isvtk) {
1433: PetscCall(DMMoabView_VTK(dm, viewer));
1434: }
1435: PetscFunctionReturn(PETSC_SUCCESS);
1436: }
1438: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
1439: {
1440: PetscFunctionBegin;
1441: dm->ops->view = DMView_Moab;
1442: dm->ops->load = NULL /* DMLoad_Moab */;
1443: dm->ops->setfromoptions = DMSetFromOptions_Moab;
1444: dm->ops->clone = DMClone_Moab;
1445: dm->ops->setup = DMSetUp_Moab;
1446: dm->ops->createlocalsection = NULL;
1447: dm->ops->createdefaultconstraints = NULL;
1448: dm->ops->createglobalvector = DMCreateGlobalVector_Moab;
1449: dm->ops->createlocalvector = DMCreateLocalVector_Moab;
1450: dm->ops->getlocaltoglobalmapping = NULL;
1451: dm->ops->createfieldis = NULL;
1452: dm->ops->createcoordinatedm = NULL /* DMCreateCoordinateDM_Moab */;
1453: dm->ops->getcoloring = NULL;
1454: dm->ops->creatematrix = DMCreateMatrix_Moab;
1455: dm->ops->createinterpolation = DMCreateInterpolation_Moab;
1456: dm->ops->createinjection = NULL /* DMCreateInjection_Moab */;
1457: dm->ops->refine = DMRefine_Moab;
1458: dm->ops->coarsen = DMCoarsen_Moab;
1459: dm->ops->refinehierarchy = DMRefineHierarchy_Moab;
1460: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Moab;
1461: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Moab;
1462: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Moab;
1463: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Moab;
1464: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Moab;
1465: dm->ops->destroy = DMDestroy_Moab;
1466: dm->ops->createsubdm = NULL /* DMCreateSubDM_Moab */;
1467: dm->ops->getdimpoints = NULL /* DMGetDimPoints_Moab */;
1468: dm->ops->locatepoints = NULL /* DMLocatePoints_Moab */;
1469: PetscFunctionReturn(PETSC_SUCCESS);
1470: }
1472: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1473: {
1474: PetscFunctionBegin;
1475: /* get all the necessary handles from the private DM object */
1476: (*newdm)->data = (DM_Moab *)dm->data;
1477: ((DM_Moab *)dm->data)->refct++;
1479: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMMOAB));
1480: PetscCall(DMInitialize_Moab(*newdm));
1481: PetscFunctionReturn(PETSC_SUCCESS);
1482: }
1484: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1485: {
1486: PetscFunctionBegin;
1488: PetscCall(PetscNew((DM_Moab **)&dm->data));
1490: ((DM_Moab *)dm->data)->bs = 1;
1491: ((DM_Moab *)dm->data)->numFields = 1;
1492: ((DM_Moab *)dm->data)->n = 0;
1493: ((DM_Moab *)dm->data)->nloc = 0;
1494: ((DM_Moab *)dm->data)->nghost = 0;
1495: ((DM_Moab *)dm->data)->nele = 0;
1496: ((DM_Moab *)dm->data)->neleloc = 0;
1497: ((DM_Moab *)dm->data)->neleghost = 0;
1498: ((DM_Moab *)dm->data)->ltog_map = NULL;
1499: ((DM_Moab *)dm->data)->ltog_sendrecv = NULL;
1501: ((DM_Moab *)dm->data)->refct = 1;
1502: ((DM_Moab *)dm->data)->parent = NULL;
1503: ((DM_Moab *)dm->data)->vlocal = new moab::Range();
1504: ((DM_Moab *)dm->data)->vowned = new moab::Range();
1505: ((DM_Moab *)dm->data)->vghost = new moab::Range();
1506: ((DM_Moab *)dm->data)->elocal = new moab::Range();
1507: ((DM_Moab *)dm->data)->eghost = new moab::Range();
1509: PetscCall(DMInitialize_Moab(dm));
1510: PetscFunctionReturn(PETSC_SUCCESS);
1511: }