Actual source code: dmmoab.cxx
1: #include <petsc/private/dmmbimpl.h>
3: #include <petscdmmoab.h>
4: #include <MBTagConventions.hpp>
5: #include <moab/NestedRefine.hpp>
6: #include <moab/Skinner.hpp>
8: /*MC
9: DMMOAB = "moab" - A DM object that encapsulates an unstructured mesh described by the MOAB mesh database.
10: Direct access to the MOAB Interface and other mesh manipulation related objects are available
11: through public API. Ability to create global and local representation of Vecs containing all
12: unknowns in the interior and shared boundary via a transparent tag-data wrapper is provided
13: along with utility functions to traverse the mesh and assemble a discrete system via
14: field-based/blocked Vec(Get/Set) methods. Input from and output to different formats are
15: available.
17: Reference: https://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html
19: Level: intermediate
21: .seealso: `DMType`, `DMMoabCreate()`, `DMCreate()`, `DMSetType()`, `DMMoabCreateMoab()`
22: M*/
24: /* External function declarations here */
25: PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
26: PETSC_EXTERN PetscErrorCode DMCreateDefaultConstraints_Moab(DM dm);
27: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
28: PETSC_EXTERN PetscErrorCode DMCreateCoordinateDM_Moab(DM dm, DM *cdm);
29: PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM *dmRefined);
30: PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM *dmCoarsened);
31: PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmRefined[]);
32: PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmCoarsened[]);
33: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm);
34: PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM, Vec *);
35: PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM, Vec *);
36: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
37: PETSC_EXTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM, Vec, InsertMode, Vec);
38: PETSC_EXTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM, Vec, InsertMode, Vec);
39: PETSC_EXTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM, Vec, InsertMode, Vec);
40: PETSC_EXTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM, Vec, InsertMode, Vec);
42: /* Un-implemented routines */
43: /*
44: PETSC_EXTERN PetscErrorCode DMCreatelocalsection_Moab(DM dm);
45: PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dmCoarse, DM dmFine, Mat *mat);
46: PETSC_EXTERN PetscErrorCode DMLoad_Moab(DM dm, PetscViewer viewer);
47: PETSC_EXTERN PetscErrorCode DMGetDimPoints_Moab(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd);
48: PETSC_EXTERN PetscErrorCode DMCreateSubDM_Moab(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
49: PETSC_EXTERN PetscErrorCode DMLocatePoints_Moab(DM dm, Vec v, IS *cellIS);
50: */
52: /*@C
53: DMMoabCreate - Creates a DMMoab object, which encapsulates a moab instance
55: Collective
57: Input Parameter:
58: . comm - The communicator for the DMMoab object
60: Output Parameter:
61: . dmb - The DMMoab object
63: Level: beginner
65: .seealso: `DMMoabCreateMoab()`
66: @*/
67: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
68: {
69: PetscFunctionBegin;
70: PetscAssertPointer(dmb, 2);
71: PetscCall(DMCreate(comm, dmb));
72: PetscCall(DMSetType(*dmb, DMMOAB));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: /*@C
77: DMMoabCreateMoab - Creates a DMMoab object, optionally from an instance and other data
79: Collective
81: Input Parameters:
82: + comm - The communicator for the DMMoab object
83: . mbiface - (ptr to) the MOAB Instance; if passed in NULL, MOAB instance is created inside PETSc, and destroyed
84: along with the DMMoab
85: . ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
86: - range - If non-NULL, contains range of entities to which DOFs will be assigned
88: Output Parameter:
89: . dmb - The DMMoab object
91: Level: intermediate
93: .seealso: `DMMoabCreate()`
94: @*/
95: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
96: {
97: moab::ErrorCode merr;
98: DM dmmb;
99: DM_Moab *dmmoab;
101: PetscFunctionBegin;
102: PetscAssertPointer(dmb, 6);
104: PetscCall(DMMoabCreate(comm, &dmmb));
105: dmmoab = (DM_Moab *)(dmmb)->data;
107: if (!mbiface) {
108: dmmoab->mbiface = new moab::Core();
109: dmmoab->icreatedinstance = PETSC_TRUE;
110: } else {
111: dmmoab->mbiface = mbiface;
112: dmmoab->icreatedinstance = PETSC_FALSE;
113: }
115: /* by default the fileset = root set. This set stores the hierarchy of entities belonging to current DM */
116: dmmoab->fileset = 0;
117: dmmoab->hlevel = 0;
118: dmmoab->nghostrings = 0;
120: #ifdef MOAB_HAVE_MPI
121: moab::EntityHandle partnset;
123: /* Create root sets for each mesh. Then pass these
124: to the load_file functions to be populated. */
125: merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, partnset);
126: MBERR("Creating partition set failed", merr);
128: /* Create the parallel communicator object with the partition handle associated with MOAB */
129: dmmoab->pcomm = moab::ParallelComm::get_pcomm(dmmoab->mbiface, partnset, &comm);
130: #endif
132: /* do the remaining initializations for DMMoab */
133: dmmoab->bs = 1;
134: dmmoab->numFields = 1;
135: PetscCall(PetscMalloc(dmmoab->numFields * sizeof(char *), &dmmoab->fieldNames));
136: PetscCall(PetscStrallocpy("DEFAULT", (char **)&dmmoab->fieldNames[0]));
137: dmmoab->rw_dbglevel = 0;
138: dmmoab->partition_by_rank = PETSC_FALSE;
139: dmmoab->extra_read_options[0] = '\0';
140: dmmoab->extra_write_options[0] = '\0';
141: dmmoab->read_mode = READ_PART;
142: dmmoab->write_mode = WRITE_PART;
144: /* set global ID tag handle */
145: if (ltog_tag && *ltog_tag) {
146: PetscCall(DMMoabSetLocalToGlobalTag(dmmb, *ltog_tag));
147: } else {
148: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
149: MBERRNM(merr);
150: if (ltog_tag) *ltog_tag = dmmoab->ltog_tag;
151: }
153: merr = dmmoab->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dmmoab->material_tag);
154: MBERRNM(merr);
156: /* set the local range of entities (vertices) of interest */
157: if (range) PetscCall(DMMoabSetLocalVertices(dmmb, range));
158: *dmb = dmmb;
159: PetscFunctionReturn(PETSC_SUCCESS);
160: }
162: #ifdef MOAB_HAVE_MPI
164: /*@C
165: DMMoabGetParallelComm - Get the ParallelComm used with this DMMoab
167: Collective
169: Input Parameter:
170: . dm - The DMMoab object being set
172: Output Parameter:
173: . pcomm - The ParallelComm for the DMMoab
175: Level: beginner
177: .seealso: `DMMoabSetInterface()`
178: @*/
179: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
180: {
181: PetscFunctionBegin;
183: *pcomm = ((DM_Moab *)(dm)->data)->pcomm;
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: #endif /* MOAB_HAVE_MPI */
189: /*@C
190: DMMoabSetInterface - Set the MOAB instance used with this DMMoab
192: Collective
194: Input Parameters:
195: + dm - The DMMoab object being set
196: - mbiface - The MOAB instance being set on this DMMoab
198: Level: beginner
200: .seealso: `DMMoabGetInterface()`
201: @*/
202: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
203: {
204: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
206: PetscFunctionBegin;
208: PetscAssertPointer(mbiface, 2);
209: #ifdef MOAB_HAVE_MPI
210: dmmoab->pcomm = NULL;
211: #endif
212: dmmoab->mbiface = mbiface;
213: dmmoab->icreatedinstance = PETSC_FALSE;
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: /*@C
218: DMMoabGetInterface - Get the MOAB instance used with this DMMoab
220: Collective
222: Input Parameter:
223: . dm - The DMMoab object being set
225: Output Parameter:
226: . mbiface - The MOAB instance set on this DMMoab
228: Level: beginner
230: .seealso: `DMMoabSetInterface()`
231: @*/
232: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
233: {
234: static PetscBool cite = PETSC_FALSE;
236: PetscFunctionBegin;
238: 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, "
239: "K. and Stimpson, C. and Ernst, C.},\n year = {2004}, note = {Report}\n}\n",
240: &cite));
241: *mbiface = ((DM_Moab *)dm->data)->mbiface;
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: /*@C
246: DMMoabSetLocalVertices - Set the entities having DOFs on this DMMoab
248: Collective
250: Input Parameters:
251: + dm - The DMMoab object being set
252: - range - The entities treated by this DMMoab
254: Level: beginner
256: .seealso: `DMMoabGetAllVertices()`
257: @*/
258: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
259: {
260: moab::Range tmpvtxs;
261: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
263: PetscFunctionBegin;
265: dmmoab->vlocal->clear();
266: dmmoab->vowned->clear();
268: dmmoab->vlocal->insert(range->begin(), range->end());
270: #ifdef MOAB_HAVE_MPI
271: moab::ErrorCode merr;
272: /* filter based on parallel status */
273: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
274: MBERRNM(merr);
276: /* filter all the non-owned and shared entities out of the list */
277: tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
278: merr = dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
279: MBERRNM(merr);
280: tmpvtxs = moab::subtract(tmpvtxs, *dmmoab->vghost);
281: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
282: #else
283: *dmmoab->vowned = *dmmoab->vlocal;
284: #endif
286: /* compute and cache the sizes of local and ghosted entities */
287: dmmoab->nloc = dmmoab->vowned->size();
288: dmmoab->nghost = dmmoab->vghost->size();
289: #ifdef MOAB_HAVE_MPI
290: PetscCall(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
291: #else
292: dmmoab->n = dmmoab->nloc;
293: #endif
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: /*@C
298: DMMoabGetAllVertices - Get the entities having DOFs on this DMMoab
300: Collective
302: Input Parameter:
303: . dm - The DMMoab object being set
305: Output Parameter:
306: . local - The local vertex entities in this DMMoab = (owned+ghosted)
308: Level: beginner
310: .seealso: `DMMoabGetLocalVertices()`
311: @*/
312: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
313: {
314: PetscFunctionBegin;
316: if (local) *local = *((DM_Moab *)dm->data)->vlocal;
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: /*@C
321: DMMoabGetLocalVertices - Get the entities having DOFs on this DMMoab
323: Collective
325: Input Parameter:
326: . dm - The DMMoab object being set
328: Output Parameters:
329: + owned - The owned vertex entities in this DMMoab
330: - ghost - The ghosted entities (non-owned) stored locally in this partition
332: Level: beginner
334: .seealso: `DMMoabGetAllVertices()`
335: @*/
336: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
337: {
338: PetscFunctionBegin;
340: if (owned) *owned = ((DM_Moab *)dm->data)->vowned;
341: if (ghost) *ghost = ((DM_Moab *)dm->data)->vghost;
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: /*@C
346: DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned
348: Collective
350: Input Parameter:
351: . dm - The DMMoab object being set
353: Output Parameter:
354: . range - The entities owned locally
356: Level: beginner
358: .seealso: `DMMoabSetLocalElements()`
359: @*/
360: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
361: {
362: PetscFunctionBegin;
364: if (range) *range = ((DM_Moab *)dm->data)->elocal;
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: /*@C
369: DMMoabSetLocalElements - Set the entities having DOFs on this DMMoab
371: Collective
373: Input Parameters:
374: + dm - The DMMoab object being set
375: - range - The entities treated by this DMMoab
377: Level: beginner
379: .seealso: `DMMoabGetLocalElements()`
380: @*/
381: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
382: {
383: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
385: PetscFunctionBegin;
387: dmmoab->elocal->clear();
388: dmmoab->eghost->clear();
389: dmmoab->elocal->insert(range->begin(), range->end());
390: #ifdef MOAB_HAVE_MPI
391: moab::ErrorCode merr;
392: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
393: MBERRNM(merr);
394: *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
395: #endif
396: dmmoab->neleloc = dmmoab->elocal->size();
397: dmmoab->neleghost = dmmoab->eghost->size();
398: #ifdef MOAB_HAVE_MPI
399: PetscCall(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
400: PetscCall(PetscInfo(dm, "Created %" PetscInt_FMT " local and %" PetscInt_FMT " global elements.\n", dmmoab->neleloc, dmmoab->nele));
401: #else
402: dmmoab->nele = dmmoab->neleloc;
403: #endif
404: PetscFunctionReturn(PETSC_SUCCESS);
405: }
407: /*@C
408: DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
410: Collective
412: Input Parameters:
413: + dm - The DMMoab object being set
414: - ltogtag - The MOAB tag used for local to global ids
416: Level: beginner
418: .seealso: `DMMoabGetLocalToGlobalTag()`
419: @*/
420: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
421: {
422: PetscFunctionBegin;
424: ((DM_Moab *)dm->data)->ltog_tag = ltogtag;
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: /*@C
429: DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
431: Collective
433: Input Parameter:
434: . dm - The DMMoab object being set
436: Output Parameter:
437: . ltog_tag - The MOAB tag used for local to global ids
439: Level: beginner
441: .seealso: `DMMoabSetLocalToGlobalTag()`
442: @*/
443: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
444: {
445: PetscFunctionBegin;
447: *ltog_tag = ((DM_Moab *)dm->data)->ltog_tag;
448: PetscFunctionReturn(PETSC_SUCCESS);
449: }
451: /*@C
452: DMMoabSetBlockSize - Set the block size used with this DMMoab
454: Collective
456: Input Parameters:
457: + dm - The DMMoab object being set
458: - bs - The block size used with this DMMoab
460: Level: beginner
462: .seealso: `DMMoabGetBlockSize()`
463: @*/
464: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
465: {
466: PetscFunctionBegin;
468: ((DM_Moab *)dm->data)->bs = bs;
469: PetscFunctionReturn(PETSC_SUCCESS);
470: }
472: /*@C
473: DMMoabGetBlockSize - Get the block size used with this DMMoab
475: Collective
477: Input Parameter:
478: . dm - The DMMoab object being set
480: Output Parameter:
481: . bs - The block size used with this DMMoab
483: Level: beginner
485: .seealso: `DMMoabSetBlockSize()`
486: @*/
487: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
488: {
489: PetscFunctionBegin;
491: *bs = ((DM_Moab *)dm->data)->bs;
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: /*@C
496: DMMoabGetSize - Get the global vertex size used with this DMMoab
498: Collective on dm
500: Input Parameter:
501: . dm - The DMMoab object being set
503: Output Parameters:
504: + neg - The number of global elements in the DMMoab instance
505: - nvg - The number of global vertices in the DMMoab instance
507: Level: beginner
509: .seealso: `DMMoabGetLocalSize()`
510: @*/
511: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
512: {
513: PetscFunctionBegin;
515: if (neg) *neg = ((DM_Moab *)dm->data)->nele;
516: if (nvg) *nvg = ((DM_Moab *)dm->data)->n;
517: PetscFunctionReturn(PETSC_SUCCESS);
518: }
520: /*@C
521: DMMoabGetLocalSize - Get the local and ghosted vertex size used with this DMMoab
523: Collective on dm
525: Input Parameter:
526: . dm - The DMMoab object being set
528: Output Parameters:
529: + nel - The number of owned elements in this processor
530: . neg - The number of ghosted elements in this processor
531: . nvl - The number of owned vertices in this processor
532: - nvg - The number of ghosted vertices in this processor
534: Level: beginner
536: .seealso: `DMMoabGetSize()`
537: @*/
538: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
539: {
540: PetscFunctionBegin;
542: if (nel) *nel = ((DM_Moab *)dm->data)->neleloc;
543: if (neg) *neg = ((DM_Moab *)dm->data)->neleghost;
544: if (nvl) *nvl = ((DM_Moab *)dm->data)->nloc;
545: if (nvg) *nvg = ((DM_Moab *)dm->data)->nghost;
546: PetscFunctionReturn(PETSC_SUCCESS);
547: }
549: /*@C
550: DMMoabGetOffset - Get the local offset for the global vector
552: Collective
554: Input Parameter:
555: . dm - The DMMoab object being set
557: Output Parameter:
558: . offset - The local offset for the global vector
560: Level: beginner
562: .seealso: `DMMoabGetDimension()`
563: @*/
564: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
565: {
566: PetscFunctionBegin;
568: *offset = ((DM_Moab *)dm->data)->vstart;
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: /*@C
573: DMMoabGetDimension - Get the dimension of the DM Mesh
575: Collective
577: Input Parameter:
578: . dm - The DMMoab object
580: Output Parameter:
581: . dim - The dimension of DM
583: Level: beginner
585: .seealso: `DMMoabGetOffset()`
586: @*/
587: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
588: {
589: PetscFunctionBegin;
591: *dim = ((DM_Moab *)dm->data)->dim;
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }
595: /*@C
596: DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
597: generated through uniform refinement.
599: Collective on dm
601: Input Parameter:
602: . dm - The DMMoab object being set
604: Output Parameter:
605: . nlevel - The current mesh hierarchy level
607: Level: beginner
609: .seealso: `DMMoabGetMaterialBlock()`
610: @*/
611: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
612: {
613: PetscFunctionBegin;
615: if (nlevel) *nlevel = ((DM_Moab *)dm->data)->hlevel;
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: /*@C
620: DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh
622: Collective
624: Input Parameters:
625: + dm - The DMMoab object
626: - ehandle - The element entity handle
628: Output Parameter:
629: . mat - The material ID for the current entity
631: Level: beginner
633: .seealso: `DMMoabGetHierarchyLevel()`
634: @*/
635: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
636: {
637: DM_Moab *dmmoab;
639: PetscFunctionBegin;
641: if (*mat) {
642: dmmoab = (DM_Moab *)(dm)->data;
643: *mat = dmmoab->materials[dmmoab->elocal->index(ehandle)];
644: }
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
648: /*@C
649: DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities
651: Collective
653: Input Parameters:
654: + dm - The DMMoab object
655: . nconn - Number of entities whose coordinates are needed
656: - conn - The vertex entity handles
658: Output Parameter:
659: . vpos - The coordinates of the requested vertex entities
661: Level: beginner
663: .seealso: `DMMoabGetVertexConnectivity()`
664: @*/
665: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
666: {
667: DM_Moab *dmmoab;
668: moab::ErrorCode merr;
670: PetscFunctionBegin;
672: PetscAssertPointer(conn, 3);
673: PetscAssertPointer(vpos, 4);
674: dmmoab = (DM_Moab *)(dm)->data;
676: /* Get connectivity information in MOAB canonical ordering */
677: if (dmmoab->hlevel) {
678: merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle *>(conn), nconn, dmmoab->hlevel, vpos);
679: MBERRNM(merr);
680: } else {
681: merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);
682: MBERRNM(merr);
683: }
684: PetscFunctionReturn(PETSC_SUCCESS);
685: }
687: /*@C
688: DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity
690: Collective
692: Input Parameters:
693: + dm - The DMMoab object
694: - vhandle - Vertex entity handle
696: Output Parameters:
697: + nconn - Number of entities whose coordinates are needed
698: - conn - The vertex entity handles
700: Level: beginner
702: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabRestoreVertexConnectivity()`
703: @*/
704: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt *nconn, moab::EntityHandle **conn)
705: {
706: DM_Moab *dmmoab;
707: std::vector<moab::EntityHandle> adj_entities, connect;
708: moab::ErrorCode merr;
710: PetscFunctionBegin;
712: PetscAssertPointer(conn, 4);
713: dmmoab = (DM_Moab *)(dm)->data;
715: /* Get connectivity information in MOAB canonical ordering */
716: merr = dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION);
717: MBERRNM(merr);
718: merr = dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect);
719: MBERRNM(merr);
721: if (conn) {
722: PetscCall(PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn));
723: PetscCall(PetscArraycpy(*conn, &connect[0], connect.size()));
724: }
725: if (nconn) *nconn = connect.size();
726: PetscFunctionReturn(PETSC_SUCCESS);
727: }
729: /*@C
730: DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity
732: Collective
734: Input Parameters:
735: + dm - The DMMoab object
736: . ehandle - Vertex entity handle
737: . nconn - Number of entities whose coordinates are needed
738: - conn - The vertex entity handles
740: Level: beginner
742: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`
743: @*/
744: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, moab::EntityHandle **conn)
745: {
746: PetscFunctionBegin;
748: PetscAssertPointer(conn, 4);
750: if (conn) PetscCall(PetscFree(*conn));
751: if (nconn) *nconn = 0;
752: PetscFunctionReturn(PETSC_SUCCESS);
753: }
755: /*@C
756: DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity
758: Collective
760: Input Parameters:
761: + dm - The DMMoab object
762: - ehandle - Vertex entity handle
764: Output Parameters:
765: + nconn - Number of entities whose coordinates are needed
766: - conn - The vertex entity handles
768: Level: beginner
770: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`, `DMMoabRestoreVertexConnectivity()`
771: @*/
772: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, const moab::EntityHandle **conn)
773: {
774: DM_Moab *dmmoab;
775: const moab::EntityHandle *connect;
776: std::vector<moab::EntityHandle> vconn;
777: moab::ErrorCode merr;
778: PetscInt nnodes;
780: PetscFunctionBegin;
782: PetscAssertPointer(conn, 4);
783: dmmoab = (DM_Moab *)(dm)->data;
785: /* Get connectivity information in MOAB canonical ordering */
786: merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes);
787: MBERRNM(merr);
788: if (conn) *conn = connect;
789: if (nconn) *nconn = nnodes;
790: PetscFunctionReturn(PETSC_SUCCESS);
791: }
793: /*@C
794: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
796: Collective
798: Input Parameters:
799: + dm - The DMMoab object
800: - ent - Entity handle
802: Output Parameter:
803: . ent_on_boundary - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
805: Level: beginner
807: .seealso: `DMMoabCheckBoundaryVertices()`
808: @*/
809: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool *ent_on_boundary)
810: {
811: moab::EntityType etype;
812: DM_Moab *dmmoab;
813: PetscInt edim;
815: PetscFunctionBegin;
817: PetscAssertPointer(ent_on_boundary, 3);
818: dmmoab = (DM_Moab *)(dm)->data;
820: /* get the entity type and handle accordingly */
821: etype = dmmoab->mbiface->type_from_handle(ent);
822: PetscCheck(etype < moab::MBPOLYHEDRON, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Entity type on the boundary skin is invalid. EntityType = %" PetscInt_FMT, etype);
824: /* get the entity dimension */
825: edim = dmmoab->mbiface->dimension_from_handle(ent);
827: *ent_on_boundary = PETSC_FALSE;
828: if (etype == moab::MBVERTEX && edim == 0) {
829: *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
830: } else {
831: if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
832: if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
833: } else { /* next check the lower-dimensional faces */
834: if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
835: }
836: }
837: PetscFunctionReturn(PETSC_SUCCESS);
838: }
840: /*@C
841: DMMoabCheckBoundaryVertices - Check whether a given entity is on the boundary (vertex, edge, face, element)
843: Input Parameters:
844: + dm - The DMMoab object
845: . nconn - Number of handles
846: - cnt - Array of entity handles
848: Output Parameter:
849: . isbdvtx - Array of boundary markers - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
851: Level: beginner
853: .seealso: `DMMoabIsEntityOnBoundary()`
854: @*/
855: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool *isbdvtx)
856: {
857: DM_Moab *dmmoab;
858: PetscInt i;
860: PetscFunctionBegin;
862: PetscAssertPointer(cnt, 3);
863: PetscAssertPointer(isbdvtx, 4);
864: dmmoab = (DM_Moab *)(dm)->data;
866: for (i = 0; i < nconn; ++i) isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
867: PetscFunctionReturn(PETSC_SUCCESS);
868: }
870: /*@C
871: DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary
873: Input Parameter:
874: . dm - The DMMoab object
876: Output Parameters:
877: + bdvtx - Boundary vertices
878: . bdelems - Boundary elements
879: - bdfaces - Boundary faces
881: Level: beginner
883: .seealso: `DMMoabCheckBoundaryVertices()`, `DMMoabIsEntityOnBoundary()`
884: @*/
885: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range **bdelems, const moab::Range **bdfaces)
886: {
887: DM_Moab *dmmoab;
889: PetscFunctionBegin;
891: dmmoab = (DM_Moab *)(dm)->data;
893: if (bdvtx) *bdvtx = dmmoab->bndyvtx;
894: if (bdfaces) *bdfaces = dmmoab->bndyfaces;
895: if (bdelems) *bdfaces = dmmoab->bndyelems;
896: PetscFunctionReturn(PETSC_SUCCESS);
897: }
899: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
900: {
901: PetscInt i;
902: moab::ErrorCode merr;
903: DM_Moab *dmmoab = (DM_Moab *)dm->data;
905: PetscFunctionBegin;
908: dmmoab->refct--;
909: if (!dmmoab->refct) {
910: delete dmmoab->vlocal;
911: delete dmmoab->vowned;
912: delete dmmoab->vghost;
913: delete dmmoab->elocal;
914: delete dmmoab->eghost;
915: delete dmmoab->bndyvtx;
916: delete dmmoab->bndyfaces;
917: delete dmmoab->bndyelems;
919: PetscCall(PetscFree(dmmoab->gsindices));
920: PetscCall(PetscFree2(dmmoab->gidmap, dmmoab->lidmap));
921: PetscCall(PetscFree(dmmoab->dfill));
922: PetscCall(PetscFree(dmmoab->ofill));
923: PetscCall(PetscFree(dmmoab->materials));
924: if (dmmoab->fieldNames) {
925: for (i = 0; i < dmmoab->numFields; i++) PetscCall(PetscFree(dmmoab->fieldNames[i]));
926: PetscCall(PetscFree(dmmoab->fieldNames));
927: }
929: if (dmmoab->nhlevels) {
930: PetscCall(PetscFree(dmmoab->hsets));
931: dmmoab->nhlevels = 0;
932: if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
933: dmmoab->hierarchy = NULL;
934: }
936: if (dmmoab->icreatedinstance) {
937: delete dmmoab->pcomm;
938: merr = dmmoab->mbiface->delete_mesh();
939: MBERRNM(merr);
940: delete dmmoab->mbiface;
941: }
942: dmmoab->mbiface = NULL;
943: #ifdef MOAB_HAVE_MPI
944: dmmoab->pcomm = NULL;
945: #endif
946: PetscCall(VecScatterDestroy(&dmmoab->ltog_sendrecv));
947: PetscCall(ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map));
948: PetscCall(PetscFree(dm->data));
949: }
950: PetscFunctionReturn(PETSC_SUCCESS);
951: }
953: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(DM dm, PetscOptionItems *PetscOptionsObject)
954: {
955: DM_Moab *dmmoab = (DM_Moab *)dm->data;
957: PetscFunctionBegin;
958: PetscOptionsHeadBegin(PetscOptionsObject, "DMMoab Options");
959: PetscCall(PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL, 0));
960: 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));
961: /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
962: 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));
963: 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));
964: PetscCall(PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum *)&dmmoab->read_mode, NULL));
965: PetscCall(PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum *)&dmmoab->write_mode, NULL));
966: PetscOptionsHeadEnd();
967: PetscFunctionReturn(PETSC_SUCCESS);
968: }
970: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
971: {
972: moab::ErrorCode merr;
973: Vec local, global;
974: IS from, to;
975: moab::Range::iterator iter;
976: PetscInt i, j, f, bs, vent, totsize, *lgmap;
977: DM_Moab *dmmoab = (DM_Moab *)dm->data;
978: moab::Range adjs;
980: PetscFunctionBegin;
982: /* Get the local and shared vertices and cache it */
983: PetscCheck(dmmoab->mbiface != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB Interface before calling SetUp.");
984: #ifdef MOAB_HAVE_MPI
985: PetscCheck(dmmoab->pcomm != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB ParallelComm object before calling SetUp.");
986: #endif
988: /* Get the entities recursively in the current part of the mesh, if user did not set the local vertices explicitly */
989: if (dmmoab->vlocal->empty()) {
990: //merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,*dmmoab->vlocal,true);MBERRNM(merr);
991: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, 0, *dmmoab->vlocal, false);
992: MBERRNM(merr);
994: #ifdef MOAB_HAVE_MPI
995: /* filter based on parallel status */
996: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
997: MBERRNM(merr);
999: /* filter all the non-owned and shared entities out of the list */
1000: // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
1001: adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
1002: merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
1003: MBERRNM(merr);
1004: adjs = moab::subtract(adjs, *dmmoab->vghost);
1005: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
1006: #else
1007: *dmmoab->vowned = *dmmoab->vlocal;
1008: #endif
1010: /* compute and cache the sizes of local and ghosted entities */
1011: dmmoab->nloc = dmmoab->vowned->size();
1012: dmmoab->nghost = dmmoab->vghost->size();
1014: #ifdef MOAB_HAVE_MPI
1015: PetscCall(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
1016: 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));
1017: #else
1018: dmmoab->n = dmmoab->nloc;
1019: #endif
1020: }
1022: {
1023: /* get the information about the local elements in the mesh */
1024: dmmoab->eghost->clear();
1026: /* first decipher the leading dimension */
1027: for (i = 3; i > 0; i--) {
1028: dmmoab->elocal->clear();
1029: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false);
1030: MBERRNM(merr);
1032: /* store the current mesh dimension */
1033: if (dmmoab->elocal->size()) {
1034: dmmoab->dim = i;
1035: break;
1036: }
1037: }
1039: PetscCall(DMSetDimension(dm, dmmoab->dim));
1041: #ifdef MOAB_HAVE_MPI
1042: /* filter the ghosted and owned element list */
1043: *dmmoab->eghost = *dmmoab->elocal;
1044: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1045: MBERRNM(merr);
1046: *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
1047: #endif
1049: dmmoab->neleloc = dmmoab->elocal->size();
1050: dmmoab->neleghost = dmmoab->eghost->size();
1052: #ifdef MOAB_HAVE_MPI
1053: PetscCall(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
1054: PetscCall(PetscInfo(NULL, "%d-dim elements: owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost));
1055: #else
1056: dmmoab->nele = dmmoab->neleloc;
1057: #endif
1058: }
1060: bs = dmmoab->bs;
1061: if (!dmmoab->ltog_tag) {
1062: /* Get the global ID tag. The global ID tag is applied to each
1063: vertex. It acts as an global identifier which MOAB uses to
1064: assemble the individual pieces of the mesh */
1065: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
1066: MBERRNM(merr);
1067: }
1069: totsize = dmmoab->vlocal->size();
1070: 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);
1071: PetscCall(PetscCalloc1(totsize, &dmmoab->gsindices));
1072: {
1073: /* first get the local indices */
1074: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]);
1075: MBERRNM(merr);
1076: if (dmmoab->nghost) { /* next get the ghosted indices */
1077: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]);
1078: MBERRNM(merr);
1079: }
1081: /* find out the local and global minima of GLOBAL_ID */
1082: dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1083: for (i = 0; i < totsize; ++i) {
1084: if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1085: if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1086: }
1088: PetscCall(MPIU_Allreduce(&dmmoab->lminmax[0], &dmmoab->gminmax[0], 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm));
1089: PetscCall(MPIU_Allreduce(&dmmoab->lminmax[1], &dmmoab->gminmax[1], 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm));
1091: /* set the GID map */
1092: for (i = 0; i < totsize; ++i) { dmmoab->gsindices[i] -= dmmoab->gminmax[0]; /* zero based index needed for IS */ }
1093: dmmoab->lminmax[0] -= dmmoab->gminmax[0];
1094: dmmoab->lminmax[1] -= dmmoab->gminmax[0];
1096: 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]));
1097: }
1098: 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,
1099: dmmoab->numFields);
1101: {
1102: dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1103: dmmoab->seqend = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1104: PetscCall(PetscInfo(NULL, "SEQUENCE: Local [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "]\n", dmmoab->seqstart, dmmoab->seqend));
1106: PetscCall(PetscMalloc2(dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->gidmap, dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->lidmap));
1107: PetscCall(PetscMalloc1(totsize * dmmoab->numFields, &lgmap));
1109: i = j = 0;
1110: /* set the owned vertex data first */
1111: for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1112: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1113: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1114: dmmoab->lidmap[vent] = i;
1115: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1116: }
1117: /* next arrange all the ghosted data information */
1118: for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1119: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1120: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1121: dmmoab->lidmap[vent] = i;
1122: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1123: }
1125: /* We need to create the Global to Local Vector Scatter Contexts
1126: 1) First create a local and global vector
1127: 2) Create a local and global IS
1128: 3) Create VecScatter and LtoGMapping objects
1129: 4) Cleanup the IS and Vec objects
1130: */
1131: PetscCall(DMCreateGlobalVector(dm, &global));
1132: PetscCall(DMCreateLocalVector(dm, &local));
1134: PetscCall(VecGetOwnershipRange(global, &dmmoab->vstart, &dmmoab->vend));
1136: /* global to local must retrieve ghost points */
1137: PetscCall(ISCreateStride(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, dmmoab->vstart, 1, &from));
1138: PetscCall(ISSetBlockSize(from, bs));
1140: PetscCall(ISCreateGeneral(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, &lgmap[0], PETSC_COPY_VALUES, &to));
1141: PetscCall(ISSetBlockSize(to, bs));
1143: if (!dmmoab->ltog_map) {
1144: /* create to the local to global mapping for vectors in order to use VecSetValuesLocal */
1145: PetscCall(ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm, dmmoab->bs, totsize * dmmoab->numFields, lgmap, PETSC_COPY_VALUES, &dmmoab->ltog_map));
1146: }
1148: /* now create the scatter object from local to global vector */
1149: PetscCall(VecScatterCreate(local, from, global, to, &dmmoab->ltog_sendrecv));
1151: /* clean up IS, Vec */
1152: PetscCall(PetscFree(lgmap));
1153: PetscCall(ISDestroy(&from));
1154: PetscCall(ISDestroy(&to));
1155: PetscCall(VecDestroy(&local));
1156: PetscCall(VecDestroy(&global));
1157: }
1159: dmmoab->bndyvtx = new moab::Range();
1160: dmmoab->bndyfaces = new moab::Range();
1161: dmmoab->bndyelems = new moab::Range();
1162: /* skin the boundary and store nodes */
1163: if (!dmmoab->hlevel) {
1164: /* get the skin vertices of boundary faces for the current partition and then filter
1165: the local, boundary faces, vertices and elements alone via PSTATUS flags;
1166: this should not give us any ghosted boundary, but if user needs such a functionality
1167: it would be easy to add it based on the find_skin query below */
1168: moab::Skinner skinner(dmmoab->mbiface);
1170: /* get the entities on the skin - only the faces */
1171: merr = skinner.find_skin(dmmoab->fileset, *dmmoab->elocal, false, *dmmoab->bndyfaces, NULL, true, true, false);
1172: MBERRNM(merr); // 'false' param indicates we want faces back, not vertices
1174: #ifdef MOAB_HAVE_MPI
1175: /* filter all the non-owned and shared entities out of the list */
1176: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1177: MBERRNM(merr);
1178: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT);
1179: MBERRNM(merr);
1180: #endif
1182: /* get all the nodes via connectivity and the parent elements via adjacency information */
1183: merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false);
1184: MBERRNM(merr);
1185: merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION);
1186: MBERRNM(merr);
1187: } else {
1188: /* Let us query the hierarchy manager and get the results directly for this level */
1189: for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1190: moab::EntityHandle elemHandle = *iter;
1191: if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1192: dmmoab->bndyelems->insert(elemHandle);
1193: /* For this boundary element, query the vertices and add them to the list */
1194: std::vector<moab::EntityHandle> connect;
1195: merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect);
1196: MBERRNM(merr);
1197: for (unsigned iv = 0; iv < connect.size(); ++iv)
1198: if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv])) dmmoab->bndyvtx->insert(connect[iv]);
1199: /* Next, let us query the boundary faces and add them also to the list */
1200: std::vector<moab::EntityHandle> faces;
1201: merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim - 1, faces);
1202: MBERRNM(merr);
1203: for (unsigned ifa = 0; ifa < faces.size(); ++ifa)
1204: if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa])) dmmoab->bndyfaces->insert(faces[ifa]);
1205: }
1206: }
1207: #ifdef MOAB_HAVE_MPI
1208: /* filter all the non-owned and shared entities out of the list */
1209: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1210: MBERRNM(merr);
1211: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1212: MBERRNM(merr);
1213: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1214: MBERRNM(merr);
1215: #endif
1216: }
1217: PetscCall(PetscInfo(NULL, "Found %zu boundary vertices, %zu boundary faces and %zu boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size()));
1219: /* Get the material sets and populate the data for all locally owned elements */
1220: {
1221: PetscCall(PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials));
1222: /* Get the count of entities of particular type from dmmoab->elocal
1223: -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1224: moab::Range msets;
1225: merr = dmmoab->mbiface->get_entities_by_type_and_tag(dmmoab->fileset, moab::MBENTITYSET, &dmmoab->material_tag, NULL, 1, msets, moab::Interface::UNION);
1226: MBERRNM(merr);
1227: if (msets.size() == 0) PetscCall(PetscInfo(NULL, "No material sets found in the fileset.\n"));
1229: for (unsigned i = 0; i < msets.size(); ++i) {
1230: moab::Range msetelems;
1231: merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);
1232: MBERRNM(merr);
1233: #ifdef MOAB_HAVE_MPI
1234: /* filter all the non-owned and shared entities out of the list */
1235: merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1236: MBERRNM(merr);
1237: #endif
1239: int partID;
1240: moab::EntityHandle mset = msets[i];
1241: merr = dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID);
1242: MBERRNM(merr);
1244: for (unsigned j = 0; j < msetelems.size(); ++j) dmmoab->materials[dmmoab->elocal->index(msetelems[j])] = partID;
1245: }
1246: }
1248: PetscFunctionReturn(PETSC_SUCCESS);
1249: }
1251: /*@C
1252: DMMoabCreateVertices - Creates and adds several vertices to the primary set represented by the DM.
1254: Collective
1256: Input Parameters:
1257: + dm - The `DM` object
1258: . coords - The connectivity of the element
1259: - nverts - The number of vertices that form the element
1261: Output Parameter:
1262: . overts - The list of vertices that were created (can be NULL)
1264: Level: beginner
1266: .seealso: `DMMoabCreateSubmesh()`, `DMMoabCreateElement()`
1267: @*/
1268: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal *coords, PetscInt nverts, moab::Range *overts)
1269: {
1270: moab::ErrorCode merr;
1271: DM_Moab *dmmoab;
1272: moab::Range verts;
1274: PetscFunctionBegin;
1276: PetscAssertPointer(coords, 2);
1278: dmmoab = (DM_Moab *)dm->data;
1280: /* Insert new points */
1281: merr = dmmoab->mbiface->create_vertices(&coords[0], nverts, verts);
1282: MBERRNM(merr);
1283: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, verts);
1284: MBERRNM(merr);
1286: if (overts) *overts = verts;
1287: PetscFunctionReturn(PETSC_SUCCESS);
1288: }
1290: /*@C
1291: DMMoabCreateElement - Adds an element of specified type to the primary set represented by the DM.
1293: Collective
1295: Input Parameters:
1296: + dm - The DM object
1297: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1298: . conn - The connectivity of the element
1299: - nverts - The number of vertices that form the element
1301: Output Parameter:
1302: . oelem - The handle to the element created and added to the DM object
1304: Level: beginner
1306: .seealso: `DMMoabCreateSubmesh()`, `DMMoabCreateVertices()`
1307: @*/
1308: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle *conn, PetscInt nverts, moab::EntityHandle *oelem)
1309: {
1310: moab::ErrorCode merr;
1311: DM_Moab *dmmoab;
1312: moab::EntityHandle elem;
1314: PetscFunctionBegin;
1316: PetscAssertPointer(conn, 3);
1318: dmmoab = (DM_Moab *)dm->data;
1320: /* Insert new element */
1321: merr = dmmoab->mbiface->create_element(type, conn, nverts, elem);
1322: MBERRNM(merr);
1323: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1);
1324: MBERRNM(merr);
1326: if (oelem) *oelem = elem;
1327: PetscFunctionReturn(PETSC_SUCCESS);
1328: }
1330: /*@C
1331: DMMoabCreateSubmesh - Creates a sub-DM object with a set that contains all vertices/elements of the parent
1332: in addition to providing support for dynamic mesh modifications. This is useful for AMR calculations to
1333: create a DM object on a refined level.
1335: Collective
1337: Input Parameters:
1338: . dm - The DM object
1340: Output Parameter:
1341: . newdm - The sub DM object with updated set information
1343: Level: advanced
1345: .seealso: `DMCreate()`, `DMMoabCreateVertices()`, `DMMoabCreateElement()`
1346: @*/
1347: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1348: {
1349: DM_Moab *dmmoab;
1350: DM_Moab *ndmmoab;
1351: moab::ErrorCode merr;
1353: PetscFunctionBegin;
1356: dmmoab = (DM_Moab *)dm->data;
1358: /* Create the basic DMMoab object and keep the default parameters created by DM impls */
1359: PetscCall(DMMoabCreateMoab(((PetscObject)dm)->comm, dmmoab->mbiface, &dmmoab->ltog_tag, NULL, newdm));
1361: /* get all the necessary handles from the private DM object */
1362: ndmmoab = (DM_Moab *)(*newdm)->data;
1364: /* set the sub-mesh's parent DM reference */
1365: ndmmoab->parent = &dm;
1367: /* create a file set to associate all entities in current mesh */
1368: merr = ndmmoab->mbiface->create_meshset(moab::MESHSET_SET, ndmmoab->fileset);
1369: MBERR("Creating file set failed", merr);
1371: /* create a meshset and then add old fileset as child */
1372: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal);
1373: MBERR("Adding child vertices to parent failed", merr);
1374: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal);
1375: MBERR("Adding child elements to parent failed", merr);
1377: /* preserve the field association between the parent and sub-mesh objects */
1378: PetscCall(DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames));
1379: PetscFunctionReturn(PETSC_SUCCESS);
1380: }
1382: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1383: {
1384: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
1385: const char *name;
1386: MPI_Comm comm;
1387: PetscMPIInt size;
1389: PetscFunctionBegin;
1390: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1391: PetscCallMPI(MPI_Comm_size(comm, &size));
1392: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1393: PetscCall(PetscViewerASCIIPushTab(viewer));
1394: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimensions:\n", name, dmmoab->dim));
1395: else PetscCall(PetscViewerASCIIPrintf(viewer, "Mesh in %" PetscInt_FMT " dimensions:\n", dmmoab->dim));
1396: /* print details about the global mesh */
1397: {
1398: PetscCall(PetscViewerASCIIPushTab(viewer));
1399: PetscCall(PetscViewerASCIIPrintf(viewer, "Sizes: cells=%" PetscInt_FMT ", vertices=%" PetscInt_FMT ", blocks=%" PetscInt_FMT "\n", dmmoab->nele, dmmoab->n, dmmoab->bs));
1400: /* print boundary data */
1401: PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary trace:\n"));
1402: {
1403: PetscCall(PetscViewerASCIIPushTab(viewer));
1404: PetscCall(PetscViewerASCIIPrintf(viewer, "cells=%zu, faces=%zu, vertices=%zu\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size()));
1405: PetscCall(PetscViewerASCIIPopTab(viewer));
1406: }
1407: /* print field data */
1408: PetscCall(PetscViewerASCIIPrintf(viewer, "Fields: %" PetscInt_FMT " components\n", dmmoab->numFields));
1409: {
1410: PetscCall(PetscViewerASCIIPushTab(viewer));
1411: for (int i = 0; i < dmmoab->numFields; ++i) PetscCall(PetscViewerASCIIPrintf(viewer, "[%" PetscInt_FMT "] - %s\n", i, dmmoab->fieldNames[i]));
1412: PetscCall(PetscViewerASCIIPopTab(viewer));
1413: }
1414: PetscCall(PetscViewerASCIIPopTab(viewer));
1415: }
1416: PetscCall(PetscViewerASCIIPopTab(viewer));
1417: PetscCall(PetscViewerFlush(viewer));
1418: PetscFunctionReturn(PETSC_SUCCESS);
1419: }
1421: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1422: {
1423: PetscFunctionReturn(PETSC_SUCCESS);
1424: }
1426: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1427: {
1428: PetscFunctionReturn(PETSC_SUCCESS);
1429: }
1431: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1432: {
1433: PetscBool iascii, ishdf5, isvtk;
1435: PetscFunctionBegin;
1438: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1439: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1440: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1441: if (iascii) {
1442: PetscCall(DMMoabView_Ascii(dm, viewer));
1443: } else if (ishdf5) {
1444: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1445: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ));
1446: PetscCall(DMMoabView_HDF5(dm, viewer));
1447: PetscCall(PetscViewerPopFormat(viewer));
1448: #else
1449: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1450: #endif
1451: } else if (isvtk) {
1452: PetscCall(DMMoabView_VTK(dm, viewer));
1453: }
1454: PetscFunctionReturn(PETSC_SUCCESS);
1455: }
1457: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
1458: {
1459: PetscFunctionBegin;
1460: dm->ops->view = DMView_Moab;
1461: dm->ops->load = NULL /* DMLoad_Moab */;
1462: dm->ops->setfromoptions = DMSetFromOptions_Moab;
1463: dm->ops->clone = DMClone_Moab;
1464: dm->ops->setup = DMSetUp_Moab;
1465: dm->ops->createlocalsection = NULL;
1466: dm->ops->createdefaultconstraints = NULL;
1467: dm->ops->createglobalvector = DMCreateGlobalVector_Moab;
1468: dm->ops->createlocalvector = DMCreateLocalVector_Moab;
1469: dm->ops->getlocaltoglobalmapping = NULL;
1470: dm->ops->createfieldis = NULL;
1471: dm->ops->createcoordinatedm = NULL /* DMCreateCoordinateDM_Moab */;
1472: dm->ops->getcoloring = NULL;
1473: dm->ops->creatematrix = DMCreateMatrix_Moab;
1474: dm->ops->createinterpolation = DMCreateInterpolation_Moab;
1475: dm->ops->createinjection = NULL /* DMCreateInjection_Moab */;
1476: dm->ops->refine = DMRefine_Moab;
1477: dm->ops->coarsen = DMCoarsen_Moab;
1478: dm->ops->refinehierarchy = DMRefineHierarchy_Moab;
1479: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Moab;
1480: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Moab;
1481: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Moab;
1482: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Moab;
1483: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Moab;
1484: dm->ops->destroy = DMDestroy_Moab;
1485: dm->ops->createsubdm = NULL /* DMCreateSubDM_Moab */;
1486: dm->ops->getdimpoints = NULL /* DMGetDimPoints_Moab */;
1487: dm->ops->locatepoints = NULL /* DMLocatePoints_Moab */;
1488: PetscFunctionReturn(PETSC_SUCCESS);
1489: }
1491: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1492: {
1493: PetscFunctionBegin;
1494: /* get all the necessary handles from the private DM object */
1495: (*newdm)->data = (DM_Moab *)dm->data;
1496: ((DM_Moab *)dm->data)->refct++;
1498: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMMOAB));
1499: PetscCall(DMInitialize_Moab(*newdm));
1500: PetscFunctionReturn(PETSC_SUCCESS);
1501: }
1503: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1504: {
1505: PetscFunctionBegin;
1507: PetscCall(PetscNew((DM_Moab **)&dm->data));
1509: ((DM_Moab *)dm->data)->bs = 1;
1510: ((DM_Moab *)dm->data)->numFields = 1;
1511: ((DM_Moab *)dm->data)->n = 0;
1512: ((DM_Moab *)dm->data)->nloc = 0;
1513: ((DM_Moab *)dm->data)->nghost = 0;
1514: ((DM_Moab *)dm->data)->nele = 0;
1515: ((DM_Moab *)dm->data)->neleloc = 0;
1516: ((DM_Moab *)dm->data)->neleghost = 0;
1517: ((DM_Moab *)dm->data)->ltog_map = NULL;
1518: ((DM_Moab *)dm->data)->ltog_sendrecv = NULL;
1520: ((DM_Moab *)dm->data)->refct = 1;
1521: ((DM_Moab *)dm->data)->parent = NULL;
1522: ((DM_Moab *)dm->data)->vlocal = new moab::Range();
1523: ((DM_Moab *)dm->data)->vowned = new moab::Range();
1524: ((DM_Moab *)dm->data)->vghost = new moab::Range();
1525: ((DM_Moab *)dm->data)->elocal = new moab::Range();
1526: ((DM_Moab *)dm->data)->eghost = new moab::Range();
1528: PetscCall(DMInitialize_Moab(dm));
1529: PetscFunctionReturn(PETSC_SUCCESS);
1530: }