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 DMCreateDefaultConstraints_Moab(DM);
25: PETSC_INTERN PetscErrorCode DMCreateMatrix_Moab(DM, Mat *);
26: PETSC_INTERN PetscErrorCode DMCreateCoordinateDM_Moab(DM, DM *);
27: PETSC_INTERN PetscErrorCode DMRefine_Moab(DM, MPI_Comm, DM *);
28: PETSC_INTERN PetscErrorCode DMCoarsen_Moab(DM, MPI_Comm, DM *);
29: PETSC_INTERN PetscErrorCode DMRefineHierarchy_Moab(DM, PetscInt, DM[]);
30: PETSC_INTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM, PetscInt, DM[]);
31: PETSC_INTERN PetscErrorCode DMCreateGlobalVector_Moab(DM, Vec *);
32: PETSC_INTERN PetscErrorCode DMCreateLocalVector_Moab(DM, Vec *);
33: PETSC_INTERN PetscErrorCode DMCreateMatrix_Moab(DM, Mat *J);
34: PETSC_INTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM, Vec, InsertMode, Vec);
35: PETSC_INTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM, Vec, InsertMode, Vec);
36: PETSC_INTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM, Vec, InsertMode, Vec);
37: PETSC_INTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM, Vec, InsertMode, Vec);
39: /* Un-implemented routines */
40: /*
41: PETSC_INTERN PetscErrorCode DMCreatelocalsection_Moab(DM);
42: PETSC_INTERN PetscErrorCode DMCreateInjection_Moab(DM, DM, Mat *);
43: PETSC_INTERN PetscErrorCode DMLoad_Moab(DM, PetscViewer);
44: PETSC_INTERN PetscErrorCode DMGetDimPoints_Moab(DM, PetscInt, PetscInt *, PetscInt *);
45: PETSC_INTERN PetscErrorCode DMCreateSubDM_Moab(DM, PetscInt, PetscInt [], IS *, DM *);
46: PETSC_INTERN PetscErrorCode DMLocatePoints_Moab(DM, Vec, IS *);
47: */
49: /*@C
50: DMMoabCreate - Creates a `DMMOAB` object, which encapsulates a moab instance
52: Collective
54: Input Parameter:
55: . comm - The communicator for the `DMMOAB` object
57: Output Parameter:
58: . dmb - The `DMMOAB` object
60: Level: beginner
62: .seealso: `DMMOAB`, `DMMoabCreateMoab()`
63: @*/
64: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
65: {
66: PetscFunctionBegin;
67: PetscAssertPointer(dmb, 2);
68: PetscCall(DMCreate(comm, dmb));
69: PetscCall(DMSetType(*dmb, DMMOAB));
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: /*@C
74: DMMoabCreateMoab - Creates a `DMMOAB` object, optionally from an instance and other data
76: Collective
78: Input Parameters:
79: + comm - The communicator for the `DMMOAB` object
80: . mbiface - (ptr to) the `DMMOAB` Instance; if passed in `NULL`, MOAB instance is created inside PETSc, and destroyed
81: along with the `DMMOAB`
82: . ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
83: - range - If non-`NULL`, contains range of entities to which DOFs will be assigned
85: Output Parameter:
86: . dmb - The `DMMOAB` object
88: Level: intermediate
90: .seealso: `DMMOAB`, `DMMoabCreate()`
91: @*/
92: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
93: {
94: DM dmmb;
95: DM_Moab *dmmoab;
97: PetscFunctionBegin;
98: PetscAssertPointer(dmb, 6);
100: PetscCall(DMMoabCreate(comm, &dmmb));
101: dmmoab = (DM_Moab *)dmmb->data;
103: if (!mbiface) {
104: dmmoab->mbiface = new moab::Core();
105: dmmoab->icreatedinstance = PETSC_TRUE;
106: } else {
107: dmmoab->mbiface = mbiface;
108: dmmoab->icreatedinstance = PETSC_FALSE;
109: }
111: /* by default the fileset = root set. This set stores the hierarchy of entities belonging to current DM */
112: dmmoab->fileset = 0;
113: dmmoab->hlevel = 0;
114: dmmoab->nghostrings = 0;
116: #ifdef MOAB_HAVE_MPI
117: moab::EntityHandle partnset;
119: /* Create root sets for each mesh. Then pass these
120: to the load_file functions to be populated. */
121: PetscCallMOAB(dmmoab->mbiface->create_meshset(moab::MESHSET_SET, partnset));
123: /* Create the parallel communicator object with the partition handle associated with MOAB */
124: dmmoab->pcomm = moab::ParallelComm::get_pcomm(dmmoab->mbiface, partnset, &comm);
125: #endif
127: /* do the remaining initializations for DMMoab */
128: dmmoab->bs = 1;
129: dmmoab->numFields = 1;
130: PetscCall(PetscMalloc(dmmoab->numFields * sizeof(char *), &dmmoab->fieldNames));
131: PetscCall(PetscStrallocpy("DEFAULT", (char **)&dmmoab->fieldNames[0]));
132: dmmoab->rw_dbglevel = 0;
133: dmmoab->partition_by_rank = PETSC_FALSE;
134: dmmoab->extra_read_options[0] = '\0';
135: dmmoab->extra_write_options[0] = '\0';
136: dmmoab->read_mode = READ_PART;
137: dmmoab->write_mode = WRITE_PART;
139: /* set global ID tag handle */
140: if (ltog_tag && *ltog_tag) PetscCall(DMMoabSetLocalToGlobalTag(dmmb, *ltog_tag));
141: else {
142: PetscCallMOAB(dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag));
143: if (ltog_tag) *ltog_tag = dmmoab->ltog_tag;
144: }
146: PetscCallMOAB(dmmoab->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dmmoab->material_tag));
148: /* set the local range of entities (vertices) of interest */
149: if (range) PetscCall(DMMoabSetLocalVertices(dmmb, range));
150: *dmb = dmmb;
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: #ifdef MOAB_HAVE_MPI
156: /*@C
157: DMMoabGetParallelComm - Get the ParallelComm used with this `DMMOAB`
159: Collective
161: Input Parameter:
162: . dm - The `DMMOAB` object being set
164: Output Parameter:
165: . pcomm - The ParallelComm for the `DMMOAB`
167: Level: beginner
169: .seealso: `DMMOAB`, `DMMoabSetInterface()`
170: @*/
171: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
172: {
173: PetscFunctionBegin;
175: *pcomm = ((DM_Moab *)dm->data)->pcomm;
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: #endif /* MOAB_HAVE_MPI */
181: /*@C
182: DMMoabSetInterface - Set the MOAB instance used with this `DMMOAB`
184: Collective
186: Input Parameters:
187: + dm - The `DMMOAB` object being set
188: - mbiface - The MOAB instance being set on this `DMMOAB`
190: Level: beginner
192: .seealso: `DMMOAB`, `DMMoabGetInterface()`
193: @*/
194: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
195: {
196: DM_Moab *dmmoab = (DM_Moab *)dm->data;
198: PetscFunctionBegin;
200: PetscAssertPointer(mbiface, 2);
201: #ifdef MOAB_HAVE_MPI
202: dmmoab->pcomm = NULL;
203: #endif
204: dmmoab->mbiface = mbiface;
205: dmmoab->icreatedinstance = PETSC_FALSE;
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: /*@C
210: DMMoabGetInterface - Get the MOAB instance used with this `DMMOAB`
212: Collective
214: Input Parameter:
215: . dm - The `DMMOAB` object being set
217: Output Parameter:
218: . mbiface - The MOAB instance set on this `DMMOAB`
220: Level: beginner
222: .seealso: `DMMOAB`, `DMMoabSetInterface()`
223: @*/
224: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
225: {
226: static PetscBool cite = PETSC_FALSE;
228: PetscFunctionBegin;
230: 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, "
231: "K. and Stimpson, C. and Ernst, C.},\n year = {2004}, note = {Report}\n}\n",
232: &cite));
233: *mbiface = ((DM_Moab *)dm->data)->mbiface;
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: /*@C
238: DMMoabSetLocalVertices - Set the entities having DOFs on this `DMMOAB`
240: Collective
242: Input Parameters:
243: + dm - The `DMMOAB` object being set
244: - range - The entities treated by this `DMMOAB`
246: Level: beginner
248: .seealso: `DMMOAB`, `DMMoabGetAllVertices()`
249: @*/
250: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
251: {
252: moab::Range tmpvtxs;
253: DM_Moab *dmmoab = (DM_Moab *)dm->data;
255: PetscFunctionBegin;
257: dmmoab->vlocal->clear();
258: dmmoab->vowned->clear();
260: dmmoab->vlocal->insert(range->begin(), range->end());
262: #ifdef MOAB_HAVE_MPI
263: /* filter based on parallel status */
264: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned));
266: /* filter all the non-owned and shared entities out of the list */
267: tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
268: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost));
269: tmpvtxs = moab::subtract(tmpvtxs, *dmmoab->vghost);
270: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
271: #else
272: *dmmoab->vowned = *dmmoab->vlocal;
273: #endif
275: /* compute and cache the sizes of local and ghosted entities */
276: dmmoab->nloc = dmmoab->vowned->size();
277: dmmoab->nghost = dmmoab->vghost->size();
278: #ifdef MOAB_HAVE_MPI
279: PetscCallMPI(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
280: #else
281: dmmoab->n = dmmoab->nloc;
282: #endif
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: /*@C
287: DMMoabGetAllVertices - Get the entities having DOFs on this `DMMOAB`
289: Collective
291: Input Parameter:
292: . dm - The `DMMOAB` object being set
294: Output Parameter:
295: . local - The local vertex entities in this `DMMOAB` = (owned+ghosted)
297: Level: beginner
299: .seealso: `DMMOAB`, `DMMoabGetLocalVertices()`
300: @*/
301: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
302: {
303: PetscFunctionBegin;
305: if (local) *local = *((DM_Moab *)dm->data)->vlocal;
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: /*@C
310: DMMoabGetLocalVertices - Get the entities having DOFs on this `DMMOAB`
312: Collective
314: Input Parameter:
315: . dm - The `DMMOAB` object being set
317: Output Parameters:
318: + owned - The owned vertex entities in this `DMMOAB`
319: - ghost - The ghosted entities (non-owned) stored locally in this partition
321: Level: beginner
323: .seealso: `DMMOAB`, `DMMoabGetAllVertices()`
324: @*/
325: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
326: {
327: PetscFunctionBegin;
329: if (owned) *owned = ((DM_Moab *)dm->data)->vowned;
330: if (ghost) *ghost = ((DM_Moab *)dm->data)->vghost;
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
334: /*@C
335: DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned
337: Collective
339: Input Parameter:
340: . dm - The `DMMOAB` object being set
342: Output Parameter:
343: . range - The entities owned locally
345: Level: beginner
347: .seealso: `DMMOAB`, `DMMoabSetLocalElements()`
348: @*/
349: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
350: {
351: PetscFunctionBegin;
353: if (range) *range = ((DM_Moab *)dm->data)->elocal;
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: /*@C
358: DMMoabSetLocalElements - Set the entities having DOFs on this `DMMOAB`
360: Collective
362: Input Parameters:
363: + dm - The `DMMOAB` object being set
364: - range - The entities treated by this `DMMOAB`
366: Level: beginner
368: .seealso: `DMMOAB`, `DMMoabGetLocalElements()`
369: @*/
370: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
371: {
372: DM_Moab *dmmoab = (DM_Moab *)dm->data;
374: PetscFunctionBegin;
376: dmmoab->elocal->clear();
377: dmmoab->eghost->clear();
378: dmmoab->elocal->insert(range->begin(), range->end());
379: #ifdef MOAB_HAVE_MPI
380: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT));
381: *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
382: #endif
383: dmmoab->neleloc = dmmoab->elocal->size();
384: dmmoab->neleghost = dmmoab->eghost->size();
385: #ifdef MOAB_HAVE_MPI
386: PetscCallMPI(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
387: PetscCall(PetscInfo(dm, "Created %" PetscInt_FMT " local and %" PetscInt_FMT " global elements.\n", dmmoab->neleloc, dmmoab->nele));
388: #else
389: dmmoab->nele = dmmoab->neleloc;
390: #endif
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: /*@C
395: DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
397: Collective
399: Input Parameters:
400: + dm - The `DMMOAB` object being set
401: - ltogtag - The `DMMOAB` tag used for local to global ids
403: Level: beginner
405: .seealso: `DMMOAB`, `DMMoabGetLocalToGlobalTag()`
406: @*/
407: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
408: {
409: PetscFunctionBegin;
411: ((DM_Moab *)dm->data)->ltog_tag = ltogtag;
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: /*@C
416: DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
418: Collective
420: Input Parameter:
421: . dm - The `DMMOAB` object being set
423: Output Parameter:
424: . ltog_tag - The MOAB tag used for local to global ids
426: Level: beginner
428: .seealso: `DMMOAB`, `DMMoabSetLocalToGlobalTag()`
429: @*/
430: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
431: {
432: PetscFunctionBegin;
434: *ltog_tag = ((DM_Moab *)dm->data)->ltog_tag;
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: /*@C
439: DMMoabSetBlockSize - Set the block size used with this `DMMOAB`
441: Collective
443: Input Parameters:
444: + dm - The `DMMOAB` object being set
445: - bs - The block size used with this `DMMOAB`
447: Level: beginner
449: .seealso: `DMMOAB`, `DMMoabGetBlockSize()`
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: .seealso: `DMMOAB`, `DMMoabSetBlockSize()`
473: @*/
474: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
475: {
476: PetscFunctionBegin;
478: *bs = ((DM_Moab *)dm->data)->bs;
479: PetscFunctionReturn(PETSC_SUCCESS);
480: }
482: /*@C
483: DMMoabGetSize - Get the global vertex size used with this `DMMOAB`
485: Collective
487: Input Parameter:
488: . dm - The `DMMOAB` object being set
490: Output Parameters:
491: + neg - The number of global elements in the `DMMOAB` instance
492: - nvg - The number of global vertices in the `DMMOAB` instance
494: Level: beginner
496: .seealso: `DMMOAB`, `DMMoabGetLocalSize()`
497: @*/
498: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
499: {
500: PetscFunctionBegin;
502: if (neg) *neg = ((DM_Moab *)dm->data)->nele;
503: if (nvg) *nvg = ((DM_Moab *)dm->data)->n;
504: PetscFunctionReturn(PETSC_SUCCESS);
505: }
507: /*@C
508: DMMoabGetLocalSize - Get the local and ghosted vertex size used with this `DMMOAB`
510: Collective
512: Input Parameter:
513: . dm - The `DMMOAB` object being set
515: Output Parameters:
516: + nel - The number of owned elements in this processor
517: . neg - The number of ghosted elements in this processor
518: . nvl - The number of owned vertices in this processor
519: - nvg - The number of ghosted vertices in this processor
521: Level: beginner
523: .seealso: `DMMOAB`, `DMMoabGetSize()`
524: @*/
525: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
526: {
527: PetscFunctionBegin;
529: if (nel) *nel = ((DM_Moab *)dm->data)->neleloc;
530: if (neg) *neg = ((DM_Moab *)dm->data)->neleghost;
531: if (nvl) *nvl = ((DM_Moab *)dm->data)->nloc;
532: if (nvg) *nvg = ((DM_Moab *)dm->data)->nghost;
533: PetscFunctionReturn(PETSC_SUCCESS);
534: }
536: /*@C
537: DMMoabGetOffset - Get the local offset for the global vector
539: Collective
541: Input Parameter:
542: . dm - The `DMMOAB` object being set
544: Output Parameter:
545: . offset - The local offset for the global vector
547: Level: beginner
549: .seealso: `DMMOAB`, `DMMoabGetDimension()`
550: @*/
551: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
552: {
553: PetscFunctionBegin;
555: *offset = ((DM_Moab *)dm->data)->vstart;
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
559: /*@C
560: DMMoabGetDimension - Get the dimension of the `DM` Mesh
562: Collective
564: Input Parameter:
565: . dm - The `DMMOAB` object
567: Output Parameter:
568: . dim - The dimension of `DM`
570: Level: beginner
572: .seealso: `DMMOAB`, `DMMoabGetOffset()`
573: @*/
574: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
575: {
576: PetscFunctionBegin;
578: *dim = ((DM_Moab *)dm->data)->dim;
579: PetscFunctionReturn(PETSC_SUCCESS);
580: }
582: /*@C
583: DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
584: generated through uniform refinement.
586: Collective
588: Input Parameter:
589: . dm - The `DMMOAB` object being set
591: Output Parameter:
592: . nlevel - The current mesh hierarchy level
594: Level: beginner
596: .seealso: `DMMOAB`, `DMMoabGetMaterialBlock()`
597: @*/
598: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
599: {
600: PetscFunctionBegin;
602: if (nlevel) *nlevel = ((DM_Moab *)dm->data)->hlevel;
603: PetscFunctionReturn(PETSC_SUCCESS);
604: }
606: /*@C
607: DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh
609: Collective
611: Input Parameters:
612: + dm - The `DMMOAB` object
613: - ehandle - The element entity handle
615: Output Parameter:
616: . mat - The material ID for the current entity
618: Level: beginner
620: .seealso: `DMMOAB`, `DMMoabGetHierarchyLevel()`
621: @*/
622: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
623: {
624: DM_Moab *dmmoab;
626: PetscFunctionBegin;
628: if (*mat) {
629: dmmoab = (DM_Moab *)dm->data;
630: *mat = dmmoab->materials[dmmoab->elocal->index(ehandle)];
631: }
632: PetscFunctionReturn(PETSC_SUCCESS);
633: }
635: /*@C
636: DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities
638: Collective
640: Input Parameters:
641: + dm - The `DMMOAB` object
642: . nconn - Number of entities whose coordinates are needed
643: - conn - The vertex entity handles
645: Output Parameter:
646: . vpos - The coordinates of the requested vertex entities
648: Level: beginner
650: .seealso: `DMMOAB`, `DMMoabGetVertexConnectivity()`
651: @*/
652: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
653: {
654: DM_Moab *dmmoab;
656: PetscFunctionBegin;
658: PetscAssertPointer(conn, 3);
659: PetscAssertPointer(vpos, 4);
660: dmmoab = (DM_Moab *)dm->data;
662: /* Get connectivity information in MOAB canonical ordering */
663: if (dmmoab->hlevel) PetscCallMOAB(dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle *>(conn), nconn, dmmoab->hlevel, vpos));
664: else PetscCallMOAB(dmmoab->mbiface->get_coords(conn, nconn, vpos));
665: PetscFunctionReturn(PETSC_SUCCESS);
666: }
668: /*@C
669: DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity
671: Collective
673: Input Parameters:
674: + dm - The `DMMOAB` object
675: - vhandle - Vertex entity handle
677: Output Parameters:
678: + nconn - Number of entities whose coordinates are needed
679: - conn - The vertex entity handles
681: Level: beginner
683: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabRestoreVertexConnectivity()`
684: @*/
685: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt *nconn, moab::EntityHandle **conn)
686: {
687: DM_Moab *dmmoab;
688: std::vector<moab::EntityHandle> adj_entities, connect;
690: PetscFunctionBegin;
692: PetscAssertPointer(conn, 4);
693: dmmoab = (DM_Moab *)dm->data;
695: /* Get connectivity information in MOAB canonical ordering */
696: PetscCallMOAB(dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION));
697: PetscCallMOAB(dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect));
699: if (conn) {
700: PetscCall(PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn));
701: PetscCall(PetscArraycpy(*conn, &connect[0], connect.size()));
702: }
703: if (nconn) *nconn = connect.size();
704: PetscFunctionReturn(PETSC_SUCCESS);
705: }
707: /*@C
708: DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity
710: Collective
712: Input Parameters:
713: + dm - The `DMMOAB` object
714: . ehandle - Vertex entity handle
715: . nconn - Number of entities whose coordinates are needed
716: - conn - The vertex entity handles
718: Level: beginner
720: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`
721: @*/
722: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, moab::EntityHandle **conn)
723: {
724: PetscFunctionBegin;
726: PetscAssertPointer(conn, 4);
728: if (conn) PetscCall(PetscFree(*conn));
729: if (nconn) *nconn = 0;
730: PetscFunctionReturn(PETSC_SUCCESS);
731: }
733: /*@C
734: DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity
736: Collective
738: Input Parameters:
739: + dm - The `DMMOAB` object
740: - ehandle - Vertex entity handle
742: Output Parameters:
743: + nconn - Number of entities whose coordinates are needed
744: - conn - The vertex entity handles
746: Level: beginner
748: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`, `DMMoabRestoreVertexConnectivity()`
749: @*/
750: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, const moab::EntityHandle **conn)
751: {
752: DM_Moab *dmmoab;
753: const moab::EntityHandle *connect;
754: std::vector<moab::EntityHandle> vconn;
755: PetscInt nnodes;
757: PetscFunctionBegin;
759: PetscAssertPointer(conn, 4);
760: dmmoab = (DM_Moab *)dm->data;
762: /* Get connectivity information in MOAB canonical ordering */
763: PetscCallMOAB(dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes));
764: if (conn) *conn = connect;
765: if (nconn) *nconn = nnodes;
766: PetscFunctionReturn(PETSC_SUCCESS);
767: }
769: /*@C
770: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
772: Collective
774: Input Parameters:
775: + dm - The `DMMOAB` object
776: - ent - Entity handle
778: Output Parameter:
779: . ent_on_boundary - `PETSC_TRUE` if entity on boundary; `PETSC_FALSE` otherwise
781: Level: beginner
783: .seealso: `DMMOAB`, `DMMoabCheckBoundaryVertices()`
784: @*/
785: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool *ent_on_boundary)
786: {
787: moab::EntityType etype;
788: DM_Moab *dmmoab;
789: PetscInt edim;
791: PetscFunctionBegin;
793: PetscAssertPointer(ent_on_boundary, 3);
794: dmmoab = (DM_Moab *)dm->data;
796: /* get the entity type and handle accordingly */
797: etype = dmmoab->mbiface->type_from_handle(ent);
798: PetscCheck(etype < moab::MBPOLYHEDRON, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Entity type on the boundary skin is invalid. EntityType = %" PetscInt_FMT, etype);
800: /* get the entity dimension */
801: edim = dmmoab->mbiface->dimension_from_handle(ent);
803: *ent_on_boundary = PETSC_FALSE;
804: if (etype == moab::MBVERTEX && edim == 0) {
805: *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
806: } else {
807: if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
808: if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
809: } else { /* next check the lower-dimensional faces */
810: if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
811: }
812: }
813: PetscFunctionReturn(PETSC_SUCCESS);
814: }
816: /*@C
817: DMMoabCheckBoundaryVertices - Check whether a given entity is on the boundary (vertex, edge, face, element)
819: Input Parameters:
820: + dm - The `DMMOAB` object
821: . nconn - Number of handles
822: - cnt - Array of entity handles
824: Output Parameter:
825: . isbdvtx - Array of boundary markers - `PETSC_TRUE` if entity on boundary; `PETSC_FALSE` otherwise
827: Level: beginner
829: .seealso: `DMMOAB`, `DMMoabIsEntityOnBoundary()`
830: @*/
831: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool *isbdvtx)
832: {
833: DM_Moab *dmmoab;
834: PetscInt i;
836: PetscFunctionBegin;
838: PetscAssertPointer(cnt, 3);
839: PetscAssertPointer(isbdvtx, 4);
840: dmmoab = (DM_Moab *)dm->data;
842: for (i = 0; i < nconn; ++i) isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
843: PetscFunctionReturn(PETSC_SUCCESS);
844: }
846: /*@C
847: DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary
849: Input Parameter:
850: . dm - The `DMMOAB` object
852: Output Parameters:
853: + bdvtx - Boundary vertices
854: . bdelems - Boundary elements
855: - bdfaces - Boundary faces
857: Level: beginner
859: .seealso: `DMMOAB`, `DMMoabCheckBoundaryVertices()`, `DMMoabIsEntityOnBoundary()`
860: @*/
861: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range **bdelems, const moab::Range **bdfaces)
862: {
863: DM_Moab *dmmoab;
865: PetscFunctionBegin;
867: dmmoab = (DM_Moab *)dm->data;
869: if (bdvtx) *bdvtx = dmmoab->bndyvtx;
870: if (bdfaces) *bdfaces = dmmoab->bndyfaces;
871: if (bdelems) *bdfaces = dmmoab->bndyelems;
872: PetscFunctionReturn(PETSC_SUCCESS);
873: }
875: static PetscErrorCode DMDestroy_Moab(DM dm)
876: {
877: PetscInt i;
878: DM_Moab *dmmoab = (DM_Moab *)dm->data;
880: PetscFunctionBegin;
883: dmmoab->refct--;
884: if (!dmmoab->refct) {
885: delete dmmoab->vlocal;
886: delete dmmoab->vowned;
887: delete dmmoab->vghost;
888: delete dmmoab->elocal;
889: delete dmmoab->eghost;
890: delete dmmoab->bndyvtx;
891: delete dmmoab->bndyfaces;
892: delete dmmoab->bndyelems;
894: PetscCall(PetscFree(dmmoab->gsindices));
895: PetscCall(PetscFree2(dmmoab->gidmap, dmmoab->lidmap));
896: PetscCall(PetscFree(dmmoab->dfill));
897: PetscCall(PetscFree(dmmoab->ofill));
898: PetscCall(PetscFree(dmmoab->materials));
899: if (dmmoab->fieldNames) {
900: for (i = 0; i < dmmoab->numFields; i++) PetscCall(PetscFree(dmmoab->fieldNames[i]));
901: PetscCall(PetscFree(dmmoab->fieldNames));
902: }
904: if (dmmoab->nhlevels) {
905: PetscCall(PetscFree(dmmoab->hsets));
906: dmmoab->nhlevels = 0;
907: if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
908: dmmoab->hierarchy = NULL;
909: }
911: if (dmmoab->icreatedinstance) {
912: delete dmmoab->pcomm;
913: PetscCallMOAB(dmmoab->mbiface->delete_mesh());
914: delete dmmoab->mbiface;
915: }
916: dmmoab->mbiface = NULL;
917: #ifdef MOAB_HAVE_MPI
918: dmmoab->pcomm = NULL;
919: #endif
920: PetscCall(VecScatterDestroy(&dmmoab->ltog_sendrecv));
921: PetscCall(ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map));
922: PetscCall(PetscFree(dm->data));
923: }
924: PetscFunctionReturn(PETSC_SUCCESS);
925: }
927: static PetscErrorCode DMSetFromOptions_Moab(DM dm, PetscOptionItems PetscOptionsObject)
928: {
929: DM_Moab *dmmoab = (DM_Moab *)dm->data;
931: PetscFunctionBegin;
932: PetscOptionsHeadBegin(PetscOptionsObject, "DMMoab Options");
933: PetscCall(PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL, 0));
934: 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));
935: /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
936: 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));
937: 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));
938: PetscCall(PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum *)&dmmoab->read_mode, NULL));
939: PetscCall(PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum *)&dmmoab->write_mode, NULL));
940: PetscOptionsHeadEnd();
941: PetscFunctionReturn(PETSC_SUCCESS);
942: }
944: static PetscErrorCode DMSetUp_Moab(DM dm)
945: {
946: Vec local, global;
947: IS from, to;
948: moab::Range::iterator iter;
949: PetscInt i, j, f, bs, vent, totsize, *lgmap;
950: DM_Moab *dmmoab = (DM_Moab *)dm->data;
951: moab::Range adjs;
953: PetscFunctionBegin;
955: /* Get the local and shared vertices and cache it */
956: PetscCheck(dmmoab->mbiface != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB Interface before calling SetUp.");
957: #ifdef MOAB_HAVE_MPI
958: PetscCheck(dmmoab->pcomm != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB ParallelComm object before calling SetUp.");
959: #endif
961: /* Get the entities recursively in the current part of the mesh, if user did not set the local vertices explicitly */
962: if (dmmoab->vlocal->empty()) {
963: PetscCallMOAB(dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, 0, *dmmoab->vlocal, false));
965: #ifdef MOAB_HAVE_MPI
966: /* filter based on parallel status */
967: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned));
969: /* filter all the non-owned and shared entities out of the list */
970: // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
971: adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
972: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost));
973: adjs = moab::subtract(adjs, *dmmoab->vghost);
974: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
975: #else
976: *dmmoab->vowned = *dmmoab->vlocal;
977: #endif
979: /* compute and cache the sizes of local and ghosted entities */
980: dmmoab->nloc = dmmoab->vowned->size();
981: dmmoab->nghost = dmmoab->vghost->size();
983: #ifdef MOAB_HAVE_MPI
984: PetscCallMPI(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
985: 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));
986: #else
987: dmmoab->n = dmmoab->nloc;
988: #endif
989: }
991: {
992: /* get the information about the local elements in the mesh */
993: dmmoab->eghost->clear();
995: /* first decipher the leading dimension */
996: for (i = 3; i > 0; i--) {
997: dmmoab->elocal->clear();
998: PetscCallMOAB(dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false));
1000: /* store the current mesh dimension */
1001: if (dmmoab->elocal->size()) {
1002: dmmoab->dim = i;
1003: break;
1004: }
1005: }
1007: PetscCall(DMSetDimension(dm, dmmoab->dim));
1009: #ifdef MOAB_HAVE_MPI
1010: /* filter the ghosted and owned element list */
1011: *dmmoab->eghost = *dmmoab->elocal;
1012: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT));
1013: *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
1014: #endif
1016: dmmoab->neleloc = dmmoab->elocal->size();
1017: dmmoab->neleghost = dmmoab->eghost->size();
1019: #ifdef MOAB_HAVE_MPI
1020: PetscCallMPI(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
1021: PetscCall(PetscInfo(NULL, "%d-dim elements: owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost));
1022: #else
1023: dmmoab->nele = dmmoab->neleloc;
1024: #endif
1025: }
1027: bs = dmmoab->bs;
1028: if (!dmmoab->ltog_tag) {
1029: /* Get the global ID tag. The global ID tag is applied to each
1030: vertex. It acts as an global identifier which MOAB uses to
1031: assemble the individual pieces of the mesh */
1032: PetscCallMOAB(dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag));
1033: }
1035: totsize = dmmoab->vlocal->size();
1036: 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);
1037: PetscCall(PetscCalloc1(totsize, &dmmoab->gsindices));
1038: {
1039: /* first get the local indices */
1040: PetscCallMOAB(dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]));
1041: /* next get the ghosted indices */
1042: if (dmmoab->nghost) PetscCallMOAB(dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]));
1044: /* find out the local and global minima of GLOBAL_ID */
1045: dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1046: for (i = 0; i < totsize; ++i) {
1047: if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1048: if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1049: }
1051: PetscCallMPI(MPIU_Allreduce(&dmmoab->lminmax[0], &dmmoab->gminmax[0], 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm));
1052: PetscCallMPI(MPIU_Allreduce(&dmmoab->lminmax[1], &dmmoab->gminmax[1], 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm));
1054: /* set the GID map */
1055: for (i = 0; i < totsize; ++i) dmmoab->gsindices[i] -= dmmoab->gminmax[0]; /* zero based index needed for IS */
1056: dmmoab->lminmax[0] -= dmmoab->gminmax[0];
1057: dmmoab->lminmax[1] -= dmmoab->gminmax[0];
1059: 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]));
1060: }
1061: 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,
1062: dmmoab->numFields);
1064: {
1065: dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1066: dmmoab->seqend = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1067: PetscCall(PetscInfo(NULL, "SEQUENCE: Local [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "]\n", dmmoab->seqstart, dmmoab->seqend));
1069: PetscCall(PetscMalloc2(dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->gidmap, dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->lidmap));
1070: PetscCall(PetscMalloc1(totsize * dmmoab->numFields, &lgmap));
1072: i = j = 0;
1073: /* set the owned vertex data first */
1074: for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1075: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1076: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1077: dmmoab->lidmap[vent] = i;
1078: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1079: }
1080: /* next arrange all the ghosted data information */
1081: for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1082: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1083: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1084: dmmoab->lidmap[vent] = i;
1085: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1086: }
1088: /* We need to create the Global to Local Vector Scatter Contexts
1089: 1) First create a local and global vector
1090: 2) Create a local and global IS
1091: 3) Create VecScatter and LtoGMapping objects
1092: 4) Cleanup the IS and Vec objects
1093: */
1094: PetscCall(DMCreateGlobalVector(dm, &global));
1095: PetscCall(DMCreateLocalVector(dm, &local));
1097: PetscCall(VecGetOwnershipRange(global, &dmmoab->vstart, &dmmoab->vend));
1099: /* global to local must retrieve ghost points */
1100: PetscCall(ISCreateStride(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, dmmoab->vstart, 1, &from));
1101: PetscCall(ISSetBlockSize(from, bs));
1103: PetscCall(ISCreateGeneral(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, &lgmap[0], PETSC_COPY_VALUES, &to));
1104: PetscCall(ISSetBlockSize(to, bs));
1106: if (!dmmoab->ltog_map) {
1107: /* create to the local to global mapping for vectors in order to use VecSetValuesLocal */
1108: PetscCall(ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm, dmmoab->bs, totsize * dmmoab->numFields, lgmap, PETSC_COPY_VALUES, &dmmoab->ltog_map));
1109: }
1111: /* now create the scatter object from local to global vector */
1112: PetscCall(VecScatterCreate(local, from, global, to, &dmmoab->ltog_sendrecv));
1114: /* clean up IS, Vec */
1115: PetscCall(PetscFree(lgmap));
1116: PetscCall(ISDestroy(&from));
1117: PetscCall(ISDestroy(&to));
1118: PetscCall(VecDestroy(&local));
1119: PetscCall(VecDestroy(&global));
1120: }
1122: dmmoab->bndyvtx = new moab::Range();
1123: dmmoab->bndyfaces = new moab::Range();
1124: dmmoab->bndyelems = new moab::Range();
1125: /* skin the boundary and store nodes */
1126: if (!dmmoab->hlevel) {
1127: /* get the skin vertices of boundary faces for the current partition and then filter
1128: the local, boundary faces, vertices and elements alone via PSTATUS flags;
1129: this should not give us any ghosted boundary, but if user needs such a functionality
1130: it would be easy to add it based on the find_skin query below */
1131: moab::Skinner skinner(dmmoab->mbiface);
1133: /* get the entities on the skin - only the faces */
1134: PetscCallMOAB(skinner.find_skin(dmmoab->fileset, *dmmoab->elocal, false, *dmmoab->bndyfaces, NULL, true, true, false));
1136: #ifdef MOAB_HAVE_MPI
1137: /* filter all the non-owned and shared entities out of the list */
1138: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT));
1139: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT));
1140: #endif
1142: /* get all the nodes via connectivity and the parent elements via adjacency information */
1143: PetscCallMOAB(dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false));
1144: PetscCallMOAB(dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION));
1145: } else {
1146: /* Let us query the hierarchy manager and get the results directly for this level */
1147: for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1148: moab::EntityHandle elemHandle = *iter;
1149: if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1150: dmmoab->bndyelems->insert(elemHandle);
1151: /* For this boundary element, query the vertices and add them to the list */
1152: std::vector<moab::EntityHandle> connect;
1153: PetscCallMOAB(dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect));
1154: for (unsigned iv = 0; iv < connect.size(); ++iv)
1155: if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv])) dmmoab->bndyvtx->insert(connect[iv]);
1156: /* Next, let us query the boundary faces and add them also to the list */
1157: std::vector<moab::EntityHandle> faces;
1158: PetscCallMOAB(dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim - 1, faces));
1159: for (unsigned ifa = 0; ifa < faces.size(); ++ifa)
1160: if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa])) dmmoab->bndyfaces->insert(faces[ifa]);
1161: }
1162: }
1163: #ifdef MOAB_HAVE_MPI
1164: /* filter all the non-owned and shared entities out of the list */
1165: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT));
1166: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT));
1167: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT));
1168: #endif
1169: }
1170: PetscCall(PetscInfo(NULL, "Found %zu boundary vertices, %zu boundary faces and %zu boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size()));
1172: /* Get the material sets and populate the data for all locally owned elements */
1173: {
1174: PetscCall(PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials));
1175: /* Get the count of entities of particular type from dmmoab->elocal
1176: -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1177: moab::Range msets;
1178: PetscCallMOAB(dmmoab->mbiface->get_entities_by_type_and_tag(dmmoab->fileset, moab::MBENTITYSET, &dmmoab->material_tag, NULL, 1, msets, moab::Interface::UNION));
1179: if (msets.size() == 0) PetscCall(PetscInfo(NULL, "No material sets found in the fileset.\n"));
1181: for (unsigned i = 0; i < msets.size(); ++i) {
1182: moab::Range msetelems;
1183: PetscCallMOAB(dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true));
1184: #ifdef MOAB_HAVE_MPI
1185: /* filter all the non-owned and shared entities out of the list */
1186: PetscCallMOAB(dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT));
1187: #endif
1189: int partID;
1190: moab::EntityHandle mset = msets[i];
1191: PetscCallMOAB(dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID));
1193: for (unsigned j = 0; j < msetelems.size(); ++j) dmmoab->materials[dmmoab->elocal->index(msetelems[j])] = partID;
1194: }
1195: }
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: /*@C
1200: DMMoabCreateVertices - Creates and adds several vertices to the primary set represented by the DM.
1202: Collective
1204: Input Parameters:
1205: + dm - The `DM` object
1206: . coords - The connectivity of the element
1207: - nverts - The number of vertices that form the element
1209: Output Parameter:
1210: . overts - The list of vertices that were created (can be `NULL`)
1212: Level: beginner
1214: .seealso: `DMMOAB`, `DMMoabCreateSubmesh()`, `DMMoabCreateElement()`
1215: @*/
1216: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal *coords, PetscInt nverts, moab::Range *overts)
1217: {
1218: DM_Moab *dmmoab;
1219: moab::Range verts;
1221: PetscFunctionBegin;
1223: PetscAssertPointer(coords, 2);
1225: dmmoab = (DM_Moab *)dm->data;
1227: /* Insert new points */
1228: PetscCallMOAB(dmmoab->mbiface->create_vertices(&coords[0], nverts, verts));
1229: PetscCallMOAB(dmmoab->mbiface->add_entities(dmmoab->fileset, verts));
1230: if (overts) *overts = verts;
1231: PetscFunctionReturn(PETSC_SUCCESS);
1232: }
1234: /*@C
1235: DMMoabCreateElement - Adds an element of specified type to the primary set represented by the DM.
1237: Collective
1239: Input Parameters:
1240: + dm - The DM object
1241: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1242: . conn - The connectivity of the element
1243: - nverts - The number of vertices that form the element
1245: Output Parameter:
1246: . oelem - The handle to the element created and added to the `DM` object
1248: Level: beginner
1250: .seealso: `DMMOAB`, `DMMoabCreateSubmesh()`, `DMMoabCreateVertices()`
1251: @*/
1252: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle *conn, PetscInt nverts, moab::EntityHandle *oelem)
1253: {
1254: DM_Moab *dmmoab;
1255: moab::EntityHandle elem;
1257: PetscFunctionBegin;
1259: PetscAssertPointer(conn, 3);
1261: dmmoab = (DM_Moab *)dm->data;
1263: /* Insert new element */
1264: PetscCallMOAB(dmmoab->mbiface->create_element(type, conn, nverts, elem));
1265: PetscCallMOAB(dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1));
1266: if (oelem) *oelem = elem;
1267: PetscFunctionReturn(PETSC_SUCCESS);
1268: }
1270: /*@C
1271: DMMoabCreateSubmesh - Creates a sub-`DM` object with a set that contains all vertices/elements of the parent
1272: in addition to providing support for dynamic mesh modifications. This is useful for AMR calculations to
1273: create a DM object on a refined level.
1275: Collective
1277: Input Parameters:
1278: . dm - The `DM` object
1280: Output Parameter:
1281: . newdm - The sub `DM` object with updated set information
1283: Level: advanced
1285: .seealso: `DMMOAB`, `DMCreate()`, `DMMoabCreateVertices()`, `DMMoabCreateElement()`
1286: @*/
1287: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1288: {
1289: DM_Moab *dmmoab;
1290: DM_Moab *ndmmoab;
1292: PetscFunctionBegin;
1295: dmmoab = (DM_Moab *)dm->data;
1297: /* Create the basic DMMOAB object and keep the default parameters created by DM impls */
1298: PetscCall(DMMoabCreateMoab(((PetscObject)dm)->comm, dmmoab->mbiface, &dmmoab->ltog_tag, NULL, newdm));
1300: /* get all the necessary handles from the private DM object */
1301: ndmmoab = (DM_Moab *)(*newdm)->data;
1303: /* set the sub-mesh's parent DM reference */
1304: ndmmoab->parent = &dm;
1306: /* create a file set to associate all entities in current mesh */
1307: PetscCallMOAB(ndmmoab->mbiface->create_meshset(moab::MESHSET_SET, ndmmoab->fileset));
1309: /* create a meshset and then add old fileset as child */
1310: PetscCallMOAB(ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal));
1311: PetscCallMOAB(ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal));
1313: /* preserve the field association between the parent and sub-mesh objects */
1314: PetscCall(DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames));
1315: PetscFunctionReturn(PETSC_SUCCESS);
1316: }
1318: static PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1319: {
1320: DM_Moab *dmmoab = (DM_Moab *)dm->data;
1321: const char *name;
1322: MPI_Comm comm;
1323: PetscMPIInt size;
1325: PetscFunctionBegin;
1326: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1327: PetscCallMPI(MPI_Comm_size(comm, &size));
1328: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1329: PetscCall(PetscViewerASCIIPushTab(viewer));
1330: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimensions:\n", name, dmmoab->dim));
1331: else PetscCall(PetscViewerASCIIPrintf(viewer, "Mesh in %" PetscInt_FMT " dimensions:\n", dmmoab->dim));
1332: /* print details about the global mesh */
1333: {
1334: PetscCall(PetscViewerASCIIPushTab(viewer));
1335: PetscCall(PetscViewerASCIIPrintf(viewer, "Sizes: cells=%" PetscInt_FMT ", vertices=%" PetscInt_FMT ", blocks=%" PetscInt_FMT "\n", dmmoab->nele, dmmoab->n, dmmoab->bs));
1336: /* print boundary data */
1337: PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary trace:\n"));
1338: {
1339: PetscCall(PetscViewerASCIIPushTab(viewer));
1340: PetscCall(PetscViewerASCIIPrintf(viewer, "cells=%zu, faces=%zu, vertices=%zu\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size()));
1341: PetscCall(PetscViewerASCIIPopTab(viewer));
1342: }
1343: /* print field data */
1344: PetscCall(PetscViewerASCIIPrintf(viewer, "Fields: %" PetscInt_FMT " components\n", dmmoab->numFields));
1345: {
1346: PetscCall(PetscViewerASCIIPushTab(viewer));
1347: for (int i = 0; i < dmmoab->numFields; ++i) PetscCall(PetscViewerASCIIPrintf(viewer, "[%" PetscInt_FMT "] - %s\n", i, dmmoab->fieldNames[i]));
1348: PetscCall(PetscViewerASCIIPopTab(viewer));
1349: }
1350: PetscCall(PetscViewerASCIIPopTab(viewer));
1351: }
1352: PetscCall(PetscViewerASCIIPopTab(viewer));
1353: PetscCall(PetscViewerFlush(viewer));
1354: PetscFunctionReturn(PETSC_SUCCESS);
1355: }
1357: static PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1358: {
1359: PetscFunctionReturn(PETSC_SUCCESS);
1360: }
1362: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1363: static PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1364: {
1365: PetscFunctionReturn(PETSC_SUCCESS);
1366: }
1367: #endif
1369: static PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1370: {
1371: PetscBool isascii, ishdf5, isvtk;
1373: PetscFunctionBegin;
1376: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1377: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1378: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1379: if (isascii) {
1380: PetscCall(DMMoabView_Ascii(dm, viewer));
1381: } else if (ishdf5) {
1382: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1383: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ));
1384: PetscCall(DMMoabView_HDF5(dm, viewer));
1385: PetscCall(PetscViewerPopFormat(viewer));
1386: #else
1387: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1388: #endif
1389: } else if (isvtk) PetscCall(DMMoabView_VTK(dm, viewer));
1390: PetscFunctionReturn(PETSC_SUCCESS);
1391: }
1393: static PetscErrorCode DMInitialize_Moab(DM);
1395: static PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1396: {
1397: PetscFunctionBegin;
1398: /* get all the necessary handles from the private DM object */
1399: (*newdm)->data = (DM_Moab *)dm->data;
1400: ((DM_Moab *)dm->data)->refct++;
1402: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMMOAB));
1403: PetscCall(DMInitialize_Moab(*newdm));
1404: PetscFunctionReturn(PETSC_SUCCESS);
1405: }
1407: static PetscErrorCode DMInitialize_Moab(DM dm)
1408: {
1409: PetscFunctionBegin;
1410: dm->ops->view = DMView_Moab;
1411: dm->ops->load = NULL /* DMLoad_Moab */;
1412: dm->ops->setfromoptions = DMSetFromOptions_Moab;
1413: dm->ops->clone = DMClone_Moab;
1414: dm->ops->setup = DMSetUp_Moab;
1415: dm->ops->createlocalsection = NULL;
1416: dm->ops->createsectionpermutation = NULL;
1417: dm->ops->createdefaultconstraints = NULL;
1418: dm->ops->createglobalvector = DMCreateGlobalVector_Moab;
1419: dm->ops->createlocalvector = DMCreateLocalVector_Moab;
1420: dm->ops->getlocaltoglobalmapping = NULL;
1421: dm->ops->createfieldis = NULL;
1422: dm->ops->createcoordinatedm = NULL /* DMCreateCoordinateDM_Moab */;
1423: dm->ops->createcellcoordinatedm = NULL;
1424: dm->ops->getcoloring = NULL;
1425: dm->ops->creatematrix = DMCreateMatrix_Moab;
1426: dm->ops->createinterpolation = DMCreateInterpolation_Moab;
1427: dm->ops->createinjection = NULL /* DMCreateInjection_Moab */;
1428: dm->ops->refine = DMRefine_Moab;
1429: dm->ops->coarsen = DMCoarsen_Moab;
1430: dm->ops->refinehierarchy = DMRefineHierarchy_Moab;
1431: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Moab;
1432: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Moab;
1433: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Moab;
1434: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Moab;
1435: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Moab;
1436: dm->ops->destroy = DMDestroy_Moab;
1437: dm->ops->createsubdm = NULL /* DMCreateSubDM_Moab */;
1438: dm->ops->getdimpoints = NULL /* DMGetDimPoints_Moab */;
1439: dm->ops->locatepoints = NULL /* DMLocatePoints_Moab */;
1440: PetscFunctionReturn(PETSC_SUCCESS);
1441: }
1443: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1444: {
1445: PetscFunctionBegin;
1447: PetscCall(PetscNew((DM_Moab **)&dm->data));
1449: ((DM_Moab *)dm->data)->bs = 1;
1450: ((DM_Moab *)dm->data)->numFields = 1;
1451: ((DM_Moab *)dm->data)->n = 0;
1452: ((DM_Moab *)dm->data)->nloc = 0;
1453: ((DM_Moab *)dm->data)->nghost = 0;
1454: ((DM_Moab *)dm->data)->nele = 0;
1455: ((DM_Moab *)dm->data)->neleloc = 0;
1456: ((DM_Moab *)dm->data)->neleghost = 0;
1457: ((DM_Moab *)dm->data)->ltog_map = NULL;
1458: ((DM_Moab *)dm->data)->ltog_sendrecv = NULL;
1460: ((DM_Moab *)dm->data)->refct = 1;
1461: ((DM_Moab *)dm->data)->parent = NULL;
1462: ((DM_Moab *)dm->data)->vlocal = new moab::Range();
1463: ((DM_Moab *)dm->data)->vowned = new moab::Range();
1464: ((DM_Moab *)dm->data)->vghost = new moab::Range();
1465: ((DM_Moab *)dm->data)->elocal = new moab::Range();
1466: ((DM_Moab *)dm->data)->eghost = new moab::Range();
1468: PetscCall(DMInitialize_Moab(dm));
1469: PetscFunctionReturn(PETSC_SUCCESS);
1470: }