Actual source code: dmmoab.cxx
1: #include <petsc/private/dmmbimpl.h>
3: #include <petscdmmoab.h>
4: #include <MBTagConventions.hpp>
5: #include <moab/NestedRefine.hpp>
6: #include <moab/Skinner.hpp>
8: /*MC
9: DMMOAB = "moab" - A `DM` object that encapsulates an unstructured mesh described by the MOAB mesh database {cite}`moabwebsite`.
10: Direct access to the MOAB Interface and other mesh manipulation related objects are available
11: through public API. Ability to create global and local representation of Vecs containing all
12: unknowns in the interior and shared boundary via a transparent tag-data wrapper is provided
13: along with utility functions to traverse the mesh and assemble a discrete system via
14: field-based/blocked Vec(Get/Set) methods. Input from and output to different formats are
15: available.
17: Level: intermediate
19: .seealso: `DMMOAB`, `DMType`, `DMMoabCreate()`, `DMCreate()`, `DMSetType()`, `DMMoabCreateMoab()`
20: M*/
22: /* External function declarations here */
23: PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
24: PETSC_EXTERN PetscErrorCode DMCreateDefaultConstraints_Moab(DM dm);
25: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
26: PETSC_EXTERN PetscErrorCode DMCreateCoordinateDM_Moab(DM dm, DM *cdm);
27: PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM *dmRefined);
28: PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM *dmCoarsened);
29: PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmRefined[]);
30: PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmCoarsened[]);
31: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm);
32: PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM, Vec *);
33: PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM, Vec *);
34: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
35: PETSC_EXTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM, Vec, InsertMode, Vec);
36: PETSC_EXTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM, Vec, InsertMode, Vec);
37: PETSC_EXTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM, Vec, InsertMode, Vec);
38: PETSC_EXTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM, Vec, InsertMode, Vec);
40: /* Un-implemented routines */
41: /*
42: PETSC_EXTERN PetscErrorCode DMCreatelocalsection_Moab(DM dm);
43: PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dmCoarse, DM dmFine, Mat *mat);
44: PETSC_EXTERN PetscErrorCode DMLoad_Moab(DM dm, PetscViewer viewer);
45: PETSC_EXTERN PetscErrorCode DMGetDimPoints_Moab(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd);
46: PETSC_EXTERN PetscErrorCode DMCreateSubDM_Moab(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
47: PETSC_EXTERN PetscErrorCode DMLocatePoints_Moab(DM dm, Vec v, IS *cellIS);
48: */
50: /*@C
51: DMMoabCreate - Creates a `DMMOAB` object, which encapsulates a moab instance
53: Collective
55: Input Parameter:
56: . comm - The communicator for the `DMMOAB` object
58: Output Parameter:
59: . dmb - The `DMMOAB` object
61: Level: beginner
63: .seealso: `DMMOAB`, `DMMoabCreateMoab()`
64: @*/
65: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
66: {
67: PetscFunctionBegin;
68: PetscAssertPointer(dmb, 2);
69: PetscCall(DMCreate(comm, dmb));
70: PetscCall(DMSetType(*dmb, DMMOAB));
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: /*@C
75: DMMoabCreateMoab - Creates a `DMMOAB` object, optionally from an instance and other data
77: Collective
79: Input Parameters:
80: + comm - The communicator for the `DMMOAB` object
81: . mbiface - (ptr to) the `DMMOAB` Instance; if passed in `NULL`, MOAB instance is created inside PETSc, and destroyed
82: along with the `DMMOAB`
83: . ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
84: - range - If non-`NULL`, contains range of entities to which DOFs will be assigned
86: Output Parameter:
87: . dmb - The `DMMOAB` object
89: Level: intermediate
91: .seealso: `DMMOAB`, `DMMoabCreate()`
92: @*/
93: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
94: {
95: moab::ErrorCode merr;
96: DM dmmb;
97: DM_Moab *dmmoab;
99: PetscFunctionBegin;
100: PetscAssertPointer(dmb, 6);
102: PetscCall(DMMoabCreate(comm, &dmmb));
103: dmmoab = (DM_Moab *)dmmb->data;
105: if (!mbiface) {
106: dmmoab->mbiface = new moab::Core();
107: dmmoab->icreatedinstance = PETSC_TRUE;
108: } else {
109: dmmoab->mbiface = mbiface;
110: dmmoab->icreatedinstance = PETSC_FALSE;
111: }
113: /* by default the fileset = root set. This set stores the hierarchy of entities belonging to current DM */
114: dmmoab->fileset = 0;
115: dmmoab->hlevel = 0;
116: dmmoab->nghostrings = 0;
118: #ifdef MOAB_HAVE_MPI
119: moab::EntityHandle partnset;
121: /* Create root sets for each mesh. Then pass these
122: to the load_file functions to be populated. */
123: merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, partnset);
124: MBERR("Creating partition set failed", merr);
126: /* Create the parallel communicator object with the partition handle associated with MOAB */
127: dmmoab->pcomm = moab::ParallelComm::get_pcomm(dmmoab->mbiface, partnset, &comm);
128: #endif
130: /* do the remaining initializations for DMMoab */
131: dmmoab->bs = 1;
132: dmmoab->numFields = 1;
133: PetscCall(PetscMalloc(dmmoab->numFields * sizeof(char *), &dmmoab->fieldNames));
134: PetscCall(PetscStrallocpy("DEFAULT", (char **)&dmmoab->fieldNames[0]));
135: dmmoab->rw_dbglevel = 0;
136: dmmoab->partition_by_rank = PETSC_FALSE;
137: dmmoab->extra_read_options[0] = '\0';
138: dmmoab->extra_write_options[0] = '\0';
139: dmmoab->read_mode = READ_PART;
140: dmmoab->write_mode = WRITE_PART;
142: /* set global ID tag handle */
143: if (ltog_tag && *ltog_tag) {
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: .seealso: `DMMOAB`, `DMMoabSetInterface()`
176: @*/
177: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
178: {
179: PetscFunctionBegin;
181: *pcomm = ((DM_Moab *)dm->data)->pcomm;
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: #endif /* MOAB_HAVE_MPI */
187: /*@C
188: DMMoabSetInterface - Set the MOAB instance used with this `DMMOAB`
190: Collective
192: Input Parameters:
193: + dm - The `DMMOAB` object being set
194: - mbiface - The MOAB instance being set on this `DMMOAB`
196: Level: beginner
198: .seealso: `DMMOAB`, `DMMoabGetInterface()`
199: @*/
200: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
201: {
202: DM_Moab *dmmoab = (DM_Moab *)dm->data;
204: PetscFunctionBegin;
206: PetscAssertPointer(mbiface, 2);
207: #ifdef MOAB_HAVE_MPI
208: dmmoab->pcomm = NULL;
209: #endif
210: dmmoab->mbiface = mbiface;
211: dmmoab->icreatedinstance = PETSC_FALSE;
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: /*@C
216: DMMoabGetInterface - Get the MOAB instance used with this `DMMOAB`
218: Collective
220: Input Parameter:
221: . dm - The `DMMOAB` object being set
223: Output Parameter:
224: . mbiface - The MOAB instance set on this `DMMOAB`
226: Level: beginner
228: .seealso: `DMMOAB`, `DMMoabSetInterface()`
229: @*/
230: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
231: {
232: static PetscBool cite = PETSC_FALSE;
234: PetscFunctionBegin;
236: 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, "
237: "K. and Stimpson, C. and Ernst, C.},\n year = {2004}, note = {Report}\n}\n",
238: &cite));
239: *mbiface = ((DM_Moab *)dm->data)->mbiface;
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: /*@C
244: DMMoabSetLocalVertices - Set the entities having DOFs on this `DMMOAB`
246: Collective
248: Input Parameters:
249: + dm - The `DMMOAB` object being set
250: - range - The entities treated by this `DMMOAB`
252: Level: beginner
254: .seealso: `DMMOAB`, `DMMoabGetAllVertices()`
255: @*/
256: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
257: {
258: moab::Range tmpvtxs;
259: DM_Moab *dmmoab = (DM_Moab *)dm->data;
261: PetscFunctionBegin;
263: dmmoab->vlocal->clear();
264: dmmoab->vowned->clear();
266: dmmoab->vlocal->insert(range->begin(), range->end());
268: #ifdef MOAB_HAVE_MPI
269: moab::ErrorCode merr;
270: /* filter based on parallel status */
271: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
272: MBERRNM(merr);
274: /* filter all the non-owned and shared entities out of the list */
275: tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
276: merr = dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
277: MBERRNM(merr);
278: tmpvtxs = moab::subtract(tmpvtxs, *dmmoab->vghost);
279: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
280: #else
281: *dmmoab->vowned = *dmmoab->vlocal;
282: #endif
284: /* compute and cache the sizes of local and ghosted entities */
285: dmmoab->nloc = dmmoab->vowned->size();
286: dmmoab->nghost = dmmoab->vghost->size();
287: #ifdef MOAB_HAVE_MPI
288: PetscCallMPI(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
289: #else
290: dmmoab->n = dmmoab->nloc;
291: #endif
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: /*@C
296: DMMoabGetAllVertices - Get the entities having DOFs on this `DMMOAB`
298: Collective
300: Input Parameter:
301: . dm - The `DMMOAB` object being set
303: Output Parameter:
304: . local - The local vertex entities in this `DMMOAB` = (owned+ghosted)
306: Level: beginner
308: .seealso: `DMMOAB`, `DMMoabGetLocalVertices()`
309: @*/
310: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
311: {
312: PetscFunctionBegin;
314: if (local) *local = *((DM_Moab *)dm->data)->vlocal;
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: /*@C
319: DMMoabGetLocalVertices - Get the entities having DOFs on this `DMMOAB`
321: Collective
323: Input Parameter:
324: . dm - The `DMMOAB` object being set
326: Output Parameters:
327: + owned - The owned vertex entities in this `DMMOAB`
328: - ghost - The ghosted entities (non-owned) stored locally in this partition
330: Level: beginner
332: .seealso: `DMMOAB`, `DMMoabGetAllVertices()`
333: @*/
334: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
335: {
336: PetscFunctionBegin;
338: if (owned) *owned = ((DM_Moab *)dm->data)->vowned;
339: if (ghost) *ghost = ((DM_Moab *)dm->data)->vghost;
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: /*@C
344: DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned
346: Collective
348: Input Parameter:
349: . dm - The `DMMOAB` object being set
351: Output Parameter:
352: . range - The entities owned locally
354: Level: beginner
356: .seealso: `DMMOAB`, `DMMoabSetLocalElements()`
357: @*/
358: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
359: {
360: PetscFunctionBegin;
362: if (range) *range = ((DM_Moab *)dm->data)->elocal;
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }
366: /*@C
367: DMMoabSetLocalElements - Set the entities having DOFs on this `DMMOAB`
369: Collective
371: Input Parameters:
372: + dm - The `DMMOAB` object being set
373: - range - The entities treated by this `DMMOAB`
375: Level: beginner
377: .seealso: `DMMOAB`, `DMMoabGetLocalElements()`
378: @*/
379: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
380: {
381: DM_Moab *dmmoab = (DM_Moab *)dm->data;
383: PetscFunctionBegin;
385: dmmoab->elocal->clear();
386: dmmoab->eghost->clear();
387: dmmoab->elocal->insert(range->begin(), range->end());
388: #ifdef MOAB_HAVE_MPI
389: moab::ErrorCode merr;
390: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
391: MBERRNM(merr);
392: *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
393: #endif
394: dmmoab->neleloc = dmmoab->elocal->size();
395: dmmoab->neleghost = dmmoab->eghost->size();
396: #ifdef MOAB_HAVE_MPI
397: PetscCallMPI(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
398: PetscCall(PetscInfo(dm, "Created %" PetscInt_FMT " local and %" PetscInt_FMT " global elements.\n", dmmoab->neleloc, dmmoab->nele));
399: #else
400: dmmoab->nele = dmmoab->neleloc;
401: #endif
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: /*@C
406: DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
408: Collective
410: Input Parameters:
411: + dm - The `DMMOAB` object being set
412: - ltogtag - The `DMMOAB` tag used for local to global ids
414: Level: beginner
416: .seealso: `DMMOAB`, `DMMoabGetLocalToGlobalTag()`
417: @*/
418: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
419: {
420: PetscFunctionBegin;
422: ((DM_Moab *)dm->data)->ltog_tag = ltogtag;
423: PetscFunctionReturn(PETSC_SUCCESS);
424: }
426: /*@C
427: DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
429: Collective
431: Input Parameter:
432: . dm - The `DMMOAB` object being set
434: Output Parameter:
435: . ltog_tag - The MOAB tag used for local to global ids
437: Level: beginner
439: .seealso: `DMMOAB`, `DMMoabSetLocalToGlobalTag()`
440: @*/
441: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
442: {
443: PetscFunctionBegin;
445: *ltog_tag = ((DM_Moab *)dm->data)->ltog_tag;
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }
449: /*@C
450: DMMoabSetBlockSize - Set the block size used with this `DMMOAB`
452: Collective
454: Input Parameters:
455: + dm - The `DMMOAB` object being set
456: - bs - The block size used with this `DMMOAB`
458: Level: beginner
460: .seealso: `DMMOAB`, `DMMoabGetBlockSize()`
461: @*/
462: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
463: {
464: PetscFunctionBegin;
466: ((DM_Moab *)dm->data)->bs = bs;
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
470: /*@C
471: DMMoabGetBlockSize - Get the block size used with this `DMMOAB`
473: Collective
475: Input Parameter:
476: . dm - The `DMMOAB` object being set
478: Output Parameter:
479: . bs - The block size used with this `DMMOAB`
481: Level: beginner
483: .seealso: `DMMOAB`, `DMMoabSetBlockSize()`
484: @*/
485: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
486: {
487: PetscFunctionBegin;
489: *bs = ((DM_Moab *)dm->data)->bs;
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }
493: /*@C
494: DMMoabGetSize - Get the global vertex size used with this `DMMOAB`
496: Collective
498: Input Parameter:
499: . dm - The `DMMOAB` object being set
501: Output Parameters:
502: + neg - The number of global elements in the `DMMOAB` instance
503: - nvg - The number of global vertices in the `DMMOAB` instance
505: Level: beginner
507: .seealso: `DMMOAB`, `DMMoabGetLocalSize()`
508: @*/
509: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
510: {
511: PetscFunctionBegin;
513: if (neg) *neg = ((DM_Moab *)dm->data)->nele;
514: if (nvg) *nvg = ((DM_Moab *)dm->data)->n;
515: PetscFunctionReturn(PETSC_SUCCESS);
516: }
518: /*@C
519: DMMoabGetLocalSize - Get the local and ghosted vertex size used with this `DMMOAB`
521: Collective
523: Input Parameter:
524: . dm - The `DMMOAB` object being set
526: Output Parameters:
527: + nel - The number of owned elements in this processor
528: . neg - The number of ghosted elements in this processor
529: . nvl - The number of owned vertices in this processor
530: - nvg - The number of ghosted vertices in this processor
532: Level: beginner
534: .seealso: `DMMOAB`, `DMMoabGetSize()`
535: @*/
536: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
537: {
538: PetscFunctionBegin;
540: if (nel) *nel = ((DM_Moab *)dm->data)->neleloc;
541: if (neg) *neg = ((DM_Moab *)dm->data)->neleghost;
542: if (nvl) *nvl = ((DM_Moab *)dm->data)->nloc;
543: if (nvg) *nvg = ((DM_Moab *)dm->data)->nghost;
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: /*@C
548: DMMoabGetOffset - Get the local offset for the global vector
550: Collective
552: Input Parameter:
553: . dm - The `DMMOAB` object being set
555: Output Parameter:
556: . offset - The local offset for the global vector
558: Level: beginner
560: .seealso: `DMMOAB`, `DMMoabGetDimension()`
561: @*/
562: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
563: {
564: PetscFunctionBegin;
566: *offset = ((DM_Moab *)dm->data)->vstart;
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: /*@C
571: DMMoabGetDimension - Get the dimension of the `DM` Mesh
573: Collective
575: Input Parameter:
576: . dm - The `DMMOAB` object
578: Output Parameter:
579: . dim - The dimension of `DM`
581: Level: beginner
583: .seealso: `DMMOAB`, `DMMoabGetOffset()`
584: @*/
585: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
586: {
587: PetscFunctionBegin;
589: *dim = ((DM_Moab *)dm->data)->dim;
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: /*@C
594: DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
595: generated through uniform refinement.
597: Collective
599: Input Parameter:
600: . dm - The `DMMOAB` object being set
602: Output Parameter:
603: . nlevel - The current mesh hierarchy level
605: Level: beginner
607: .seealso: `DMMOAB`, `DMMoabGetMaterialBlock()`
608: @*/
609: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
610: {
611: PetscFunctionBegin;
613: if (nlevel) *nlevel = ((DM_Moab *)dm->data)->hlevel;
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: /*@C
618: DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh
620: Collective
622: Input Parameters:
623: + dm - The `DMMOAB` object
624: - ehandle - The element entity handle
626: Output Parameter:
627: . mat - The material ID for the current entity
629: Level: beginner
631: .seealso: `DMMOAB`, `DMMoabGetHierarchyLevel()`
632: @*/
633: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
634: {
635: DM_Moab *dmmoab;
637: PetscFunctionBegin;
639: if (*mat) {
640: dmmoab = (DM_Moab *)dm->data;
641: *mat = dmmoab->materials[dmmoab->elocal->index(ehandle)];
642: }
643: PetscFunctionReturn(PETSC_SUCCESS);
644: }
646: /*@C
647: DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities
649: Collective
651: Input Parameters:
652: + dm - The `DMMOAB` object
653: . nconn - Number of entities whose coordinates are needed
654: - conn - The vertex entity handles
656: Output Parameter:
657: . vpos - The coordinates of the requested vertex entities
659: Level: beginner
661: .seealso: `DMMOAB`, `DMMoabGetVertexConnectivity()`
662: @*/
663: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
664: {
665: DM_Moab *dmmoab;
666: moab::ErrorCode merr;
668: PetscFunctionBegin;
670: PetscAssertPointer(conn, 3);
671: PetscAssertPointer(vpos, 4);
672: dmmoab = (DM_Moab *)dm->data;
674: /* Get connectivity information in MOAB canonical ordering */
675: if (dmmoab->hlevel) {
676: merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle *>(conn), nconn, dmmoab->hlevel, vpos);
677: MBERRNM(merr);
678: } else {
679: merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);
680: MBERRNM(merr);
681: }
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: /*@C
686: DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity
688: Collective
690: Input Parameters:
691: + dm - The `DMMOAB` object
692: - vhandle - Vertex entity handle
694: Output Parameters:
695: + nconn - Number of entities whose coordinates are needed
696: - conn - The vertex entity handles
698: Level: beginner
700: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabRestoreVertexConnectivity()`
701: @*/
702: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt *nconn, moab::EntityHandle **conn)
703: {
704: DM_Moab *dmmoab;
705: std::vector<moab::EntityHandle> adj_entities, connect;
706: moab::ErrorCode merr;
708: PetscFunctionBegin;
710: PetscAssertPointer(conn, 4);
711: dmmoab = (DM_Moab *)dm->data;
713: /* Get connectivity information in MOAB canonical ordering */
714: merr = dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION);
715: MBERRNM(merr);
716: merr = dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect);
717: MBERRNM(merr);
719: if (conn) {
720: PetscCall(PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn));
721: PetscCall(PetscArraycpy(*conn, &connect[0], connect.size()));
722: }
723: if (nconn) *nconn = connect.size();
724: PetscFunctionReturn(PETSC_SUCCESS);
725: }
727: /*@C
728: DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity
730: Collective
732: Input Parameters:
733: + dm - The `DMMOAB` object
734: . ehandle - Vertex entity handle
735: . nconn - Number of entities whose coordinates are needed
736: - conn - The vertex entity handles
738: Level: beginner
740: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`
741: @*/
742: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, moab::EntityHandle **conn)
743: {
744: PetscFunctionBegin;
746: PetscAssertPointer(conn, 4);
748: if (conn) PetscCall(PetscFree(*conn));
749: if (nconn) *nconn = 0;
750: PetscFunctionReturn(PETSC_SUCCESS);
751: }
753: /*@C
754: DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity
756: Collective
758: Input Parameters:
759: + dm - The `DMMOAB` object
760: - ehandle - Vertex entity handle
762: Output Parameters:
763: + nconn - Number of entities whose coordinates are needed
764: - conn - The vertex entity handles
766: Level: beginner
768: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`, `DMMoabRestoreVertexConnectivity()`
769: @*/
770: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, const moab::EntityHandle **conn)
771: {
772: DM_Moab *dmmoab;
773: const moab::EntityHandle *connect;
774: std::vector<moab::EntityHandle> vconn;
775: moab::ErrorCode merr;
776: PetscInt nnodes;
778: PetscFunctionBegin;
780: PetscAssertPointer(conn, 4);
781: dmmoab = (DM_Moab *)dm->data;
783: /* Get connectivity information in MOAB canonical ordering */
784: merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes);
785: MBERRNM(merr);
786: if (conn) *conn = connect;
787: if (nconn) *nconn = nnodes;
788: PetscFunctionReturn(PETSC_SUCCESS);
789: }
791: /*@C
792: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
794: Collective
796: Input Parameters:
797: + dm - The `DMMOAB` object
798: - ent - Entity handle
800: Output Parameter:
801: . ent_on_boundary - `PETSC_TRUE` if entity on boundary; `PETSC_FALSE` otherwise
803: Level: beginner
805: .seealso: `DMMOAB`, `DMMoabCheckBoundaryVertices()`
806: @*/
807: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool *ent_on_boundary)
808: {
809: moab::EntityType etype;
810: DM_Moab *dmmoab;
811: PetscInt edim;
813: PetscFunctionBegin;
815: PetscAssertPointer(ent_on_boundary, 3);
816: dmmoab = (DM_Moab *)dm->data;
818: /* get the entity type and handle accordingly */
819: etype = dmmoab->mbiface->type_from_handle(ent);
820: PetscCheck(etype < moab::MBPOLYHEDRON, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Entity type on the boundary skin is invalid. EntityType = %" PetscInt_FMT, etype);
822: /* get the entity dimension */
823: edim = dmmoab->mbiface->dimension_from_handle(ent);
825: *ent_on_boundary = PETSC_FALSE;
826: if (etype == moab::MBVERTEX && edim == 0) {
827: *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
828: } else {
829: if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
830: if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
831: } else { /* next check the lower-dimensional faces */
832: if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
833: }
834: }
835: PetscFunctionReturn(PETSC_SUCCESS);
836: }
838: /*@C
839: DMMoabCheckBoundaryVertices - Check whether a given entity is on the boundary (vertex, edge, face, element)
841: Input Parameters:
842: + dm - The `DMMOAB` object
843: . nconn - Number of handles
844: - cnt - Array of entity handles
846: Output Parameter:
847: . isbdvtx - Array of boundary markers - `PETSC_TRUE` if entity on boundary; `PETSC_FALSE` otherwise
849: Level: beginner
851: .seealso: `DMMOAB`, `DMMoabIsEntityOnBoundary()`
852: @*/
853: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool *isbdvtx)
854: {
855: DM_Moab *dmmoab;
856: PetscInt i;
858: PetscFunctionBegin;
860: PetscAssertPointer(cnt, 3);
861: PetscAssertPointer(isbdvtx, 4);
862: dmmoab = (DM_Moab *)dm->data;
864: for (i = 0; i < nconn; ++i) isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
865: PetscFunctionReturn(PETSC_SUCCESS);
866: }
868: /*@C
869: DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary
871: Input Parameter:
872: . dm - The `DMMOAB` object
874: Output Parameters:
875: + bdvtx - Boundary vertices
876: . bdelems - Boundary elements
877: - bdfaces - Boundary faces
879: Level: beginner
881: .seealso: `DMMOAB`, `DMMoabCheckBoundaryVertices()`, `DMMoabIsEntityOnBoundary()`
882: @*/
883: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range **bdelems, const moab::Range **bdfaces)
884: {
885: DM_Moab *dmmoab;
887: PetscFunctionBegin;
889: dmmoab = (DM_Moab *)dm->data;
891: if (bdvtx) *bdvtx = dmmoab->bndyvtx;
892: if (bdfaces) *bdfaces = dmmoab->bndyfaces;
893: if (bdelems) *bdfaces = dmmoab->bndyelems;
894: PetscFunctionReturn(PETSC_SUCCESS);
895: }
897: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
898: {
899: PetscInt i;
900: moab::ErrorCode merr;
901: DM_Moab *dmmoab = (DM_Moab *)dm->data;
903: PetscFunctionBegin;
906: dmmoab->refct--;
907: if (!dmmoab->refct) {
908: delete dmmoab->vlocal;
909: delete dmmoab->vowned;
910: delete dmmoab->vghost;
911: delete dmmoab->elocal;
912: delete dmmoab->eghost;
913: delete dmmoab->bndyvtx;
914: delete dmmoab->bndyfaces;
915: delete dmmoab->bndyelems;
917: PetscCall(PetscFree(dmmoab->gsindices));
918: PetscCall(PetscFree2(dmmoab->gidmap, dmmoab->lidmap));
919: PetscCall(PetscFree(dmmoab->dfill));
920: PetscCall(PetscFree(dmmoab->ofill));
921: PetscCall(PetscFree(dmmoab->materials));
922: if (dmmoab->fieldNames) {
923: for (i = 0; i < dmmoab->numFields; i++) PetscCall(PetscFree(dmmoab->fieldNames[i]));
924: PetscCall(PetscFree(dmmoab->fieldNames));
925: }
927: if (dmmoab->nhlevels) {
928: PetscCall(PetscFree(dmmoab->hsets));
929: dmmoab->nhlevels = 0;
930: if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
931: dmmoab->hierarchy = NULL;
932: }
934: if (dmmoab->icreatedinstance) {
935: delete dmmoab->pcomm;
936: merr = dmmoab->mbiface->delete_mesh();
937: MBERRNM(merr);
938: delete dmmoab->mbiface;
939: }
940: dmmoab->mbiface = NULL;
941: #ifdef MOAB_HAVE_MPI
942: dmmoab->pcomm = NULL;
943: #endif
944: PetscCall(VecScatterDestroy(&dmmoab->ltog_sendrecv));
945: PetscCall(ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map));
946: PetscCall(PetscFree(dm->data));
947: }
948: PetscFunctionReturn(PETSC_SUCCESS);
949: }
951: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(DM dm, PetscOptionItems PetscOptionsObject)
952: {
953: DM_Moab *dmmoab = (DM_Moab *)dm->data;
955: PetscFunctionBegin;
956: PetscOptionsHeadBegin(PetscOptionsObject, "DMMoab Options");
957: PetscCall(PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL, 0));
958: 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));
959: /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
960: 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));
961: 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));
962: PetscCall(PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum *)&dmmoab->read_mode, NULL));
963: PetscCall(PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum *)&dmmoab->write_mode, NULL));
964: PetscOptionsHeadEnd();
965: PetscFunctionReturn(PETSC_SUCCESS);
966: }
968: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
969: {
970: moab::ErrorCode merr;
971: Vec local, global;
972: IS from, to;
973: moab::Range::iterator iter;
974: PetscInt i, j, f, bs, vent, totsize, *lgmap;
975: DM_Moab *dmmoab = (DM_Moab *)dm->data;
976: moab::Range adjs;
978: PetscFunctionBegin;
980: /* Get the local and shared vertices and cache it */
981: PetscCheck(dmmoab->mbiface != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB Interface before calling SetUp.");
982: #ifdef MOAB_HAVE_MPI
983: PetscCheck(dmmoab->pcomm != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB ParallelComm object before calling SetUp.");
984: #endif
986: /* Get the entities recursively in the current part of the mesh, if user did not set the local vertices explicitly */
987: if (dmmoab->vlocal->empty()) {
988: //merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,*dmmoab->vlocal,true);MBERRNM(merr);
989: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, 0, *dmmoab->vlocal, false);
990: MBERRNM(merr);
992: #ifdef MOAB_HAVE_MPI
993: /* filter based on parallel status */
994: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
995: MBERRNM(merr);
997: /* filter all the non-owned and shared entities out of the list */
998: // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
999: adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
1000: merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
1001: MBERRNM(merr);
1002: adjs = moab::subtract(adjs, *dmmoab->vghost);
1003: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
1004: #else
1005: *dmmoab->vowned = *dmmoab->vlocal;
1006: #endif
1008: /* compute and cache the sizes of local and ghosted entities */
1009: dmmoab->nloc = dmmoab->vowned->size();
1010: dmmoab->nghost = dmmoab->vghost->size();
1012: #ifdef MOAB_HAVE_MPI
1013: PetscCallMPI(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
1014: 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));
1015: #else
1016: dmmoab->n = dmmoab->nloc;
1017: #endif
1018: }
1020: {
1021: /* get the information about the local elements in the mesh */
1022: dmmoab->eghost->clear();
1024: /* first decipher the leading dimension */
1025: for (i = 3; i > 0; i--) {
1026: dmmoab->elocal->clear();
1027: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false);
1028: MBERRNM(merr);
1030: /* store the current mesh dimension */
1031: if (dmmoab->elocal->size()) {
1032: dmmoab->dim = i;
1033: break;
1034: }
1035: }
1037: PetscCall(DMSetDimension(dm, dmmoab->dim));
1039: #ifdef MOAB_HAVE_MPI
1040: /* filter the ghosted and owned element list */
1041: *dmmoab->eghost = *dmmoab->elocal;
1042: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1043: MBERRNM(merr);
1044: *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
1045: #endif
1047: dmmoab->neleloc = dmmoab->elocal->size();
1048: dmmoab->neleghost = dmmoab->eghost->size();
1050: #ifdef MOAB_HAVE_MPI
1051: PetscCallMPI(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
1052: PetscCall(PetscInfo(NULL, "%d-dim elements: owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost));
1053: #else
1054: dmmoab->nele = dmmoab->neleloc;
1055: #endif
1056: }
1058: bs = dmmoab->bs;
1059: if (!dmmoab->ltog_tag) {
1060: /* Get the global ID tag. The global ID tag is applied to each
1061: vertex. It acts as an global identifier which MOAB uses to
1062: assemble the individual pieces of the mesh */
1063: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
1064: MBERRNM(merr);
1065: }
1067: totsize = dmmoab->vlocal->size();
1068: 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);
1069: PetscCall(PetscCalloc1(totsize, &dmmoab->gsindices));
1070: {
1071: /* first get the local indices */
1072: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]);
1073: MBERRNM(merr);
1074: if (dmmoab->nghost) { /* next get the ghosted indices */
1075: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]);
1076: MBERRNM(merr);
1077: }
1079: /* find out the local and global minima of GLOBAL_ID */
1080: dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1081: for (i = 0; i < totsize; ++i) {
1082: if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1083: if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1084: }
1086: PetscCallMPI(MPIU_Allreduce(&dmmoab->lminmax[0], &dmmoab->gminmax[0], 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm));
1087: PetscCallMPI(MPIU_Allreduce(&dmmoab->lminmax[1], &dmmoab->gminmax[1], 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm));
1089: /* set the GID map */
1090: for (i = 0; i < totsize; ++i) { dmmoab->gsindices[i] -= dmmoab->gminmax[0]; /* zero based index needed for IS */ }
1091: dmmoab->lminmax[0] -= dmmoab->gminmax[0];
1092: dmmoab->lminmax[1] -= dmmoab->gminmax[0];
1094: 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]));
1095: }
1096: 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,
1097: dmmoab->numFields);
1099: {
1100: dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1101: dmmoab->seqend = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1102: PetscCall(PetscInfo(NULL, "SEQUENCE: Local [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "]\n", dmmoab->seqstart, dmmoab->seqend));
1104: PetscCall(PetscMalloc2(dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->gidmap, dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->lidmap));
1105: PetscCall(PetscMalloc1(totsize * dmmoab->numFields, &lgmap));
1107: i = j = 0;
1108: /* set the owned vertex data first */
1109: for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1110: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1111: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1112: dmmoab->lidmap[vent] = i;
1113: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1114: }
1115: /* next arrange all the ghosted data information */
1116: for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1117: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1118: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1119: dmmoab->lidmap[vent] = i;
1120: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1121: }
1123: /* We need to create the Global to Local Vector Scatter Contexts
1124: 1) First create a local and global vector
1125: 2) Create a local and global IS
1126: 3) Create VecScatter and LtoGMapping objects
1127: 4) Cleanup the IS and Vec objects
1128: */
1129: PetscCall(DMCreateGlobalVector(dm, &global));
1130: PetscCall(DMCreateLocalVector(dm, &local));
1132: PetscCall(VecGetOwnershipRange(global, &dmmoab->vstart, &dmmoab->vend));
1134: /* global to local must retrieve ghost points */
1135: PetscCall(ISCreateStride(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, dmmoab->vstart, 1, &from));
1136: PetscCall(ISSetBlockSize(from, bs));
1138: PetscCall(ISCreateGeneral(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, &lgmap[0], PETSC_COPY_VALUES, &to));
1139: PetscCall(ISSetBlockSize(to, bs));
1141: if (!dmmoab->ltog_map) {
1142: /* create to the local to global mapping for vectors in order to use VecSetValuesLocal */
1143: PetscCall(ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm, dmmoab->bs, totsize * dmmoab->numFields, lgmap, PETSC_COPY_VALUES, &dmmoab->ltog_map));
1144: }
1146: /* now create the scatter object from local to global vector */
1147: PetscCall(VecScatterCreate(local, from, global, to, &dmmoab->ltog_sendrecv));
1149: /* clean up IS, Vec */
1150: PetscCall(PetscFree(lgmap));
1151: PetscCall(ISDestroy(&from));
1152: PetscCall(ISDestroy(&to));
1153: PetscCall(VecDestroy(&local));
1154: PetscCall(VecDestroy(&global));
1155: }
1157: dmmoab->bndyvtx = new moab::Range();
1158: dmmoab->bndyfaces = new moab::Range();
1159: dmmoab->bndyelems = new moab::Range();
1160: /* skin the boundary and store nodes */
1161: if (!dmmoab->hlevel) {
1162: /* get the skin vertices of boundary faces for the current partition and then filter
1163: the local, boundary faces, vertices and elements alone via PSTATUS flags;
1164: this should not give us any ghosted boundary, but if user needs such a functionality
1165: it would be easy to add it based on the find_skin query below */
1166: moab::Skinner skinner(dmmoab->mbiface);
1168: /* get the entities on the skin - only the faces */
1169: merr = skinner.find_skin(dmmoab->fileset, *dmmoab->elocal, false, *dmmoab->bndyfaces, NULL, true, true, false);
1170: MBERRNM(merr); // 'false' param indicates we want faces back, not vertices
1172: #ifdef MOAB_HAVE_MPI
1173: /* filter all the non-owned and shared entities out of the list */
1174: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1175: MBERRNM(merr);
1176: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT);
1177: MBERRNM(merr);
1178: #endif
1180: /* get all the nodes via connectivity and the parent elements via adjacency information */
1181: merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false);
1182: MBERRNM(merr);
1183: merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION);
1184: MBERRNM(merr);
1185: } else {
1186: /* Let us query the hierarchy manager and get the results directly for this level */
1187: for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1188: moab::EntityHandle elemHandle = *iter;
1189: if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1190: dmmoab->bndyelems->insert(elemHandle);
1191: /* For this boundary element, query the vertices and add them to the list */
1192: std::vector<moab::EntityHandle> connect;
1193: merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect);
1194: MBERRNM(merr);
1195: for (unsigned iv = 0; iv < connect.size(); ++iv)
1196: if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv])) dmmoab->bndyvtx->insert(connect[iv]);
1197: /* Next, let us query the boundary faces and add them also to the list */
1198: std::vector<moab::EntityHandle> faces;
1199: merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim - 1, faces);
1200: MBERRNM(merr);
1201: for (unsigned ifa = 0; ifa < faces.size(); ++ifa)
1202: if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa])) dmmoab->bndyfaces->insert(faces[ifa]);
1203: }
1204: }
1205: #ifdef MOAB_HAVE_MPI
1206: /* filter all the non-owned and shared entities out of the list */
1207: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1208: MBERRNM(merr);
1209: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1210: MBERRNM(merr);
1211: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1212: MBERRNM(merr);
1213: #endif
1214: }
1215: PetscCall(PetscInfo(NULL, "Found %zu boundary vertices, %zu boundary faces and %zu boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size()));
1217: /* Get the material sets and populate the data for all locally owned elements */
1218: {
1219: PetscCall(PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials));
1220: /* Get the count of entities of particular type from dmmoab->elocal
1221: -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1222: moab::Range msets;
1223: merr = dmmoab->mbiface->get_entities_by_type_and_tag(dmmoab->fileset, moab::MBENTITYSET, &dmmoab->material_tag, NULL, 1, msets, moab::Interface::UNION);
1224: MBERRNM(merr);
1225: if (msets.size() == 0) PetscCall(PetscInfo(NULL, "No material sets found in the fileset.\n"));
1227: for (unsigned i = 0; i < msets.size(); ++i) {
1228: moab::Range msetelems;
1229: merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);
1230: MBERRNM(merr);
1231: #ifdef MOAB_HAVE_MPI
1232: /* filter all the non-owned and shared entities out of the list */
1233: merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1234: MBERRNM(merr);
1235: #endif
1237: int partID;
1238: moab::EntityHandle mset = msets[i];
1239: merr = dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID);
1240: MBERRNM(merr);
1242: for (unsigned j = 0; j < msetelems.size(); ++j) dmmoab->materials[dmmoab->elocal->index(msetelems[j])] = partID;
1243: }
1244: }
1245: PetscFunctionReturn(PETSC_SUCCESS);
1246: }
1248: /*@C
1249: DMMoabCreateVertices - Creates and adds several vertices to the primary set represented by the DM.
1251: Collective
1253: Input Parameters:
1254: + dm - The `DM` object
1255: . coords - The connectivity of the element
1256: - nverts - The number of vertices that form the element
1258: Output Parameter:
1259: . overts - The list of vertices that were created (can be `NULL`)
1261: Level: beginner
1263: .seealso: `DMMOAB`, `DMMoabCreateSubmesh()`, `DMMoabCreateElement()`
1264: @*/
1265: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal *coords, PetscInt nverts, moab::Range *overts)
1266: {
1267: moab::ErrorCode merr;
1268: DM_Moab *dmmoab;
1269: moab::Range verts;
1271: PetscFunctionBegin;
1273: PetscAssertPointer(coords, 2);
1275: dmmoab = (DM_Moab *)dm->data;
1277: /* Insert new points */
1278: merr = dmmoab->mbiface->create_vertices(&coords[0], nverts, verts);
1279: MBERRNM(merr);
1280: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, verts);
1281: MBERRNM(merr);
1283: if (overts) *overts = verts;
1284: PetscFunctionReturn(PETSC_SUCCESS);
1285: }
1287: /*@C
1288: DMMoabCreateElement - Adds an element of specified type to the primary set represented by the DM.
1290: Collective
1292: Input Parameters:
1293: + dm - The DM object
1294: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1295: . conn - The connectivity of the element
1296: - nverts - The number of vertices that form the element
1298: Output Parameter:
1299: . oelem - The handle to the element created and added to the `DM` object
1301: Level: beginner
1303: .seealso: `DMMOAB`, `DMMoabCreateSubmesh()`, `DMMoabCreateVertices()`
1304: @*/
1305: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle *conn, PetscInt nverts, moab::EntityHandle *oelem)
1306: {
1307: moab::ErrorCode merr;
1308: DM_Moab *dmmoab;
1309: moab::EntityHandle elem;
1311: PetscFunctionBegin;
1313: PetscAssertPointer(conn, 3);
1315: dmmoab = (DM_Moab *)dm->data;
1317: /* Insert new element */
1318: merr = dmmoab->mbiface->create_element(type, conn, nverts, elem);
1319: MBERRNM(merr);
1320: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1);
1321: MBERRNM(merr);
1323: if (oelem) *oelem = elem;
1324: PetscFunctionReturn(PETSC_SUCCESS);
1325: }
1327: /*@C
1328: DMMoabCreateSubmesh - Creates a sub-`DM` object with a set that contains all vertices/elements of the parent
1329: in addition to providing support for dynamic mesh modifications. This is useful for AMR calculations to
1330: create a DM object on a refined level.
1332: Collective
1334: Input Parameters:
1335: . dm - The `DM` object
1337: Output Parameter:
1338: . newdm - The sub `DM` object with updated set information
1340: Level: advanced
1342: .seealso: `DMMOAB`, `DMCreate()`, `DMMoabCreateVertices()`, `DMMoabCreateElement()`
1343: @*/
1344: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1345: {
1346: DM_Moab *dmmoab;
1347: DM_Moab *ndmmoab;
1348: moab::ErrorCode merr;
1350: PetscFunctionBegin;
1353: dmmoab = (DM_Moab *)dm->data;
1355: /* Create the basic DMMOAB object and keep the default parameters created by DM impls */
1356: PetscCall(DMMoabCreateMoab(((PetscObject)dm)->comm, dmmoab->mbiface, &dmmoab->ltog_tag, NULL, newdm));
1358: /* get all the necessary handles from the private DM object */
1359: ndmmoab = (DM_Moab *)(*newdm)->data;
1361: /* set the sub-mesh's parent DM reference */
1362: ndmmoab->parent = &dm;
1364: /* create a file set to associate all entities in current mesh */
1365: merr = ndmmoab->mbiface->create_meshset(moab::MESHSET_SET, ndmmoab->fileset);
1366: MBERR("Creating file set failed", merr);
1368: /* create a meshset and then add old fileset as child */
1369: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal);
1370: MBERR("Adding child vertices to parent failed", merr);
1371: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal);
1372: MBERR("Adding child elements to parent failed", merr);
1374: /* preserve the field association between the parent and sub-mesh objects */
1375: PetscCall(DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames));
1376: PetscFunctionReturn(PETSC_SUCCESS);
1377: }
1379: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1380: {
1381: DM_Moab *dmmoab = (DM_Moab *)dm->data;
1382: const char *name;
1383: MPI_Comm comm;
1384: PetscMPIInt size;
1386: PetscFunctionBegin;
1387: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1388: PetscCallMPI(MPI_Comm_size(comm, &size));
1389: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1390: PetscCall(PetscViewerASCIIPushTab(viewer));
1391: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimensions:\n", name, dmmoab->dim));
1392: else PetscCall(PetscViewerASCIIPrintf(viewer, "Mesh in %" PetscInt_FMT " dimensions:\n", dmmoab->dim));
1393: /* print details about the global mesh */
1394: {
1395: PetscCall(PetscViewerASCIIPushTab(viewer));
1396: PetscCall(PetscViewerASCIIPrintf(viewer, "Sizes: cells=%" PetscInt_FMT ", vertices=%" PetscInt_FMT ", blocks=%" PetscInt_FMT "\n", dmmoab->nele, dmmoab->n, dmmoab->bs));
1397: /* print boundary data */
1398: PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary trace:\n"));
1399: {
1400: PetscCall(PetscViewerASCIIPushTab(viewer));
1401: PetscCall(PetscViewerASCIIPrintf(viewer, "cells=%zu, faces=%zu, vertices=%zu\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size()));
1402: PetscCall(PetscViewerASCIIPopTab(viewer));
1403: }
1404: /* print field data */
1405: PetscCall(PetscViewerASCIIPrintf(viewer, "Fields: %" PetscInt_FMT " components\n", dmmoab->numFields));
1406: {
1407: PetscCall(PetscViewerASCIIPushTab(viewer));
1408: for (int i = 0; i < dmmoab->numFields; ++i) PetscCall(PetscViewerASCIIPrintf(viewer, "[%" PetscInt_FMT "] - %s\n", i, dmmoab->fieldNames[i]));
1409: PetscCall(PetscViewerASCIIPopTab(viewer));
1410: }
1411: PetscCall(PetscViewerASCIIPopTab(viewer));
1412: }
1413: PetscCall(PetscViewerASCIIPopTab(viewer));
1414: PetscCall(PetscViewerFlush(viewer));
1415: PetscFunctionReturn(PETSC_SUCCESS);
1416: }
1418: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1419: {
1420: PetscFunctionReturn(PETSC_SUCCESS);
1421: }
1423: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1424: {
1425: PetscFunctionReturn(PETSC_SUCCESS);
1426: }
1428: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1429: {
1430: PetscBool iascii, ishdf5, isvtk;
1432: PetscFunctionBegin;
1435: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1436: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1437: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1438: if (iascii) {
1439: PetscCall(DMMoabView_Ascii(dm, viewer));
1440: } else if (ishdf5) {
1441: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1442: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ));
1443: PetscCall(DMMoabView_HDF5(dm, viewer));
1444: PetscCall(PetscViewerPopFormat(viewer));
1445: #else
1446: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1447: #endif
1448: } else if (isvtk) {
1449: PetscCall(DMMoabView_VTK(dm, viewer));
1450: }
1451: PetscFunctionReturn(PETSC_SUCCESS);
1452: }
1454: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
1455: {
1456: PetscFunctionBegin;
1457: dm->ops->view = DMView_Moab;
1458: dm->ops->load = NULL /* DMLoad_Moab */;
1459: dm->ops->setfromoptions = DMSetFromOptions_Moab;
1460: dm->ops->clone = DMClone_Moab;
1461: dm->ops->setup = DMSetUp_Moab;
1462: dm->ops->createlocalsection = NULL;
1463: dm->ops->createsectionpermutation = NULL;
1464: dm->ops->createdefaultconstraints = NULL;
1465: dm->ops->createglobalvector = DMCreateGlobalVector_Moab;
1466: dm->ops->createlocalvector = DMCreateLocalVector_Moab;
1467: dm->ops->getlocaltoglobalmapping = NULL;
1468: dm->ops->createfieldis = NULL;
1469: dm->ops->createcoordinatedm = NULL /* DMCreateCoordinateDM_Moab */;
1470: dm->ops->getcoloring = NULL;
1471: dm->ops->creatematrix = DMCreateMatrix_Moab;
1472: dm->ops->createinterpolation = DMCreateInterpolation_Moab;
1473: dm->ops->createinjection = NULL /* DMCreateInjection_Moab */;
1474: dm->ops->refine = DMRefine_Moab;
1475: dm->ops->coarsen = DMCoarsen_Moab;
1476: dm->ops->refinehierarchy = DMRefineHierarchy_Moab;
1477: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Moab;
1478: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Moab;
1479: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Moab;
1480: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Moab;
1481: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Moab;
1482: dm->ops->destroy = DMDestroy_Moab;
1483: dm->ops->createsubdm = NULL /* DMCreateSubDM_Moab */;
1484: dm->ops->getdimpoints = NULL /* DMGetDimPoints_Moab */;
1485: dm->ops->locatepoints = NULL /* DMLocatePoints_Moab */;
1486: PetscFunctionReturn(PETSC_SUCCESS);
1487: }
1489: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1490: {
1491: PetscFunctionBegin;
1492: /* get all the necessary handles from the private DM object */
1493: (*newdm)->data = (DM_Moab *)dm->data;
1494: ((DM_Moab *)dm->data)->refct++;
1496: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMMOAB));
1497: PetscCall(DMInitialize_Moab(*newdm));
1498: PetscFunctionReturn(PETSC_SUCCESS);
1499: }
1501: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1502: {
1503: PetscFunctionBegin;
1505: PetscCall(PetscNew((DM_Moab **)&dm->data));
1507: ((DM_Moab *)dm->data)->bs = 1;
1508: ((DM_Moab *)dm->data)->numFields = 1;
1509: ((DM_Moab *)dm->data)->n = 0;
1510: ((DM_Moab *)dm->data)->nloc = 0;
1511: ((DM_Moab *)dm->data)->nghost = 0;
1512: ((DM_Moab *)dm->data)->nele = 0;
1513: ((DM_Moab *)dm->data)->neleloc = 0;
1514: ((DM_Moab *)dm->data)->neleghost = 0;
1515: ((DM_Moab *)dm->data)->ltog_map = NULL;
1516: ((DM_Moab *)dm->data)->ltog_sendrecv = NULL;
1518: ((DM_Moab *)dm->data)->refct = 1;
1519: ((DM_Moab *)dm->data)->parent = NULL;
1520: ((DM_Moab *)dm->data)->vlocal = new moab::Range();
1521: ((DM_Moab *)dm->data)->vowned = new moab::Range();
1522: ((DM_Moab *)dm->data)->vghost = new moab::Range();
1523: ((DM_Moab *)dm->data)->elocal = new moab::Range();
1524: ((DM_Moab *)dm->data)->eghost = new moab::Range();
1526: PetscCall(DMInitialize_Moab(dm));
1527: PetscFunctionReturn(PETSC_SUCCESS);
1528: }