Actual source code: dmmoab.cxx
1: #include <petsc/private/dmmbimpl.h>
3: #include <petscdmmoab.h>
4: #include <MBTagConventions.hpp>
5: #include <moab/NestedRefine.hpp>
6: #include <moab/Skinner.hpp>
8: /*MC
9: DMMOAB = "moab" - A DM object that encapsulates an unstructured mesh described by the MOAB mesh database.
10: Direct access to the MOAB Interface and other mesh manipulation related objects are available
11: through public API. Ability to create global and local representation of Vecs containing all
12: unknowns in the interior and shared boundary via a transparent tag-data wrapper is provided
13: along with utility functions to traverse the mesh and assemble a discrete system via
14: field-based/blocked Vec(Get/Set) methods. Input from and output to different formats are
15: available.
17: Reference: https://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html
19: Level: intermediate
21: .seealso: `DMType`, `DMMoabCreate()`, `DMCreate()`, `DMSetType()`, `DMMoabCreateMoab()`
22: M*/
24: /* External function declarations here */
25: PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
26: PETSC_EXTERN PetscErrorCode DMCreateDefaultConstraints_Moab(DM dm);
27: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
28: PETSC_EXTERN PetscErrorCode DMCreateCoordinateDM_Moab(DM dm, DM *cdm);
29: PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM *dmRefined);
30: PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM *dmCoarsened);
31: PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmRefined[]);
32: PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmCoarsened[]);
33: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm);
34: PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM, Vec *);
35: PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM, Vec *);
36: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
37: PETSC_EXTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM, Vec, InsertMode, Vec);
38: PETSC_EXTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM, Vec, InsertMode, Vec);
39: PETSC_EXTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM, Vec, InsertMode, Vec);
40: PETSC_EXTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM, Vec, InsertMode, Vec);
42: /* Un-implemented routines */
43: /*
44: PETSC_EXTERN PetscErrorCode DMCreatelocalsection_Moab(DM dm);
45: PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dmCoarse, DM dmFine, Mat *mat);
46: PETSC_EXTERN PetscErrorCode DMLoad_Moab(DM dm, PetscViewer viewer);
47: PETSC_EXTERN PetscErrorCode DMGetDimPoints_Moab(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd);
48: PETSC_EXTERN PetscErrorCode DMCreateSubDM_Moab(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
49: PETSC_EXTERN PetscErrorCode DMLocatePoints_Moab(DM dm, Vec v, IS *cellIS);
50: */
52: /*@C
53: DMMoabCreate - Creates a DMMoab object, which encapsulates a moab instance
55: Collective
57: Input Parameter:
58: . comm - The communicator for the DMMoab object
60: Output Parameter:
61: . dmb - The DMMoab object
63: Level: beginner
65: @*/
66: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
67: {
69: DMCreate(comm, dmb);
70: DMSetType(*dmb, DMMOAB);
71: return 0;
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 MOAB 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: @*/
92: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
93: {
94: moab::ErrorCode merr;
95: DM dmmb;
96: DM_Moab *dmmoab;
100: DMMoabCreate(comm, &dmmb);
101: dmmoab = (DM_Moab *)(dmmb)->data;
103: if (!mbiface) {
104: dmmoab->mbiface = new moab::Core();
105: dmmoab->icreatedinstance = PETSC_TRUE;
106: } else {
107: dmmoab->mbiface = mbiface;
108: dmmoab->icreatedinstance = PETSC_FALSE;
109: }
111: /* by default the fileset = root set. This set stores the hierarchy of entities belonging to current DM */
112: dmmoab->fileset = 0;
113: dmmoab->hlevel = 0;
114: dmmoab->nghostrings = 0;
116: #ifdef MOAB_HAVE_MPI
117: moab::EntityHandle partnset;
119: /* Create root sets for each mesh. Then pass these
120: to the load_file functions to be populated. */
121: merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, partnset);
122: MBERR("Creating partition set failed", merr);
124: /* Create the parallel communicator object with the partition handle associated with MOAB */
125: dmmoab->pcomm = moab::ParallelComm::get_pcomm(dmmoab->mbiface, partnset, &comm);
126: #endif
128: /* do the remaining initializations for DMMoab */
129: dmmoab->bs = 1;
130: dmmoab->numFields = 1;
131: PetscMalloc(dmmoab->numFields * sizeof(char *), &dmmoab->fieldNames);
132: PetscStrallocpy("DEFAULT", (char **)&dmmoab->fieldNames[0]);
133: dmmoab->rw_dbglevel = 0;
134: dmmoab->partition_by_rank = PETSC_FALSE;
135: dmmoab->extra_read_options[0] = '\0';
136: dmmoab->extra_write_options[0] = '\0';
137: dmmoab->read_mode = READ_PART;
138: dmmoab->write_mode = WRITE_PART;
140: /* set global ID tag handle */
141: if (ltog_tag && *ltog_tag) {
142: DMMoabSetLocalToGlobalTag(dmmb, *ltog_tag);
143: } else {
144: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
145: MBERRNM(merr);
146: if (ltog_tag) *ltog_tag = dmmoab->ltog_tag;
147: }
149: merr = dmmoab->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dmmoab->material_tag);
150: MBERRNM(merr);
152: /* set the local range of entities (vertices) of interest */
153: if (range) DMMoabSetLocalVertices(dmmb, range);
154: *dmb = dmmb;
155: return 0;
156: }
158: #ifdef MOAB_HAVE_MPI
160: /*@C
161: DMMoabGetParallelComm - Get the ParallelComm used with this DMMoab
163: Collective
165: Input Parameter:
166: . dm - The DMMoab object being set
168: Output Parameter:
169: . pcomm - The ParallelComm for the DMMoab
171: Level: beginner
173: @*/
174: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
175: {
177: *pcomm = ((DM_Moab *)(dm)->data)->pcomm;
178: return 0;
179: }
181: #endif /* MOAB_HAVE_MPI */
183: /*@C
184: DMMoabSetInterface - Set the MOAB instance used with this DMMoab
186: Collective
188: Input Parameters:
189: + dm - The DMMoab object being set
190: - mbiface - The MOAB instance being set on this DMMoab
192: Level: beginner
194: @*/
195: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
196: {
197: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
201: #ifdef MOAB_HAVE_MPI
202: dmmoab->pcomm = NULL;
203: #endif
204: dmmoab->mbiface = mbiface;
205: dmmoab->icreatedinstance = PETSC_FALSE;
206: return 0;
207: }
209: /*@C
210: DMMoabGetInterface - Get the MOAB instance used with this DMMoab
212: Collective
214: Input Parameter:
215: . dm - The DMMoab object being set
217: Output Parameter:
218: . mbiface - The MOAB instance set on this DMMoab
220: Level: beginner
222: @*/
223: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
224: {
225: static PetscBool cite = PETSC_FALSE;
228: PetscCall(PetscCitationsRegister("@techreport{tautges_moab:_2004,\n type = {{SAND2004-1592}},\n title = {{MOAB:} A Mesh-Oriented Database}, institution = {Sandia National Laboratories},\n author = {Tautges, T. J. and Meyers, R. and Merkley, "
229: "K. and Stimpson, C. and Ernst, C.},\n year = {2004}, note = {Report}\n}\n",
230: &cite));
231: *mbiface = ((DM_Moab *)dm->data)->mbiface;
232: return 0;
233: }
235: /*@C
236: DMMoabSetLocalVertices - Set the entities having DOFs on this DMMoab
238: Collective
240: Input Parameters:
241: + dm - The DMMoab object being set
242: - range - The entities treated by this DMMoab
244: Level: beginner
246: @*/
247: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
248: {
249: moab::Range tmpvtxs;
250: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
253: dmmoab->vlocal->clear();
254: dmmoab->vowned->clear();
256: dmmoab->vlocal->insert(range->begin(), range->end());
258: #ifdef MOAB_HAVE_MPI
259: moab::ErrorCode merr;
260: /* filter based on parallel status */
261: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
262: MBERRNM(merr);
264: /* filter all the non-owned and shared entities out of the list */
265: tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
266: merr = dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
267: MBERRNM(merr);
268: tmpvtxs = moab::subtract(tmpvtxs, *dmmoab->vghost);
269: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
270: #else
271: *dmmoab->vowned = *dmmoab->vlocal;
272: #endif
274: /* compute and cache the sizes of local and ghosted entities */
275: dmmoab->nloc = dmmoab->vowned->size();
276: dmmoab->nghost = dmmoab->vghost->size();
277: #ifdef MOAB_HAVE_MPI
278: MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
279: #else
280: dmmoab->n = dmmoab->nloc;
281: #endif
282: return 0;
283: }
285: /*@C
286: DMMoabGetAllVertices - Get the entities having DOFs on this DMMoab
288: Collective
290: Input Parameter:
291: . dm - The DMMoab object being set
293: Output Parameter:
294: . owned - The local vertex entities in this DMMoab = (owned+ghosted)
296: Level: beginner
298: @*/
299: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
300: {
302: if (local) *local = *((DM_Moab *)dm->data)->vlocal;
303: return 0;
304: }
306: /*@C
307: DMMoabGetLocalVertices - Get the entities having DOFs on this DMMoab
309: Collective
311: Input Parameter:
312: . dm - The DMMoab object being set
314: Output Parameters:
315: + owned - The owned vertex entities in this DMMoab
316: - ghost - The ghosted entities (non-owned) stored locally in this partition
318: Level: beginner
320: @*/
321: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
322: {
324: if (owned) *owned = ((DM_Moab *)dm->data)->vowned;
325: if (ghost) *ghost = ((DM_Moab *)dm->data)->vghost;
326: return 0;
327: }
329: /*@C
330: DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned
332: Collective
334: Input Parameter:
335: . dm - The DMMoab object being set
337: Output Parameter:
338: . range - The entities owned locally
340: Level: beginner
342: @*/
343: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
344: {
346: if (range) *range = ((DM_Moab *)dm->data)->elocal;
347: return 0;
348: }
350: /*@C
351: DMMoabSetLocalElements - Set the entities having DOFs on this DMMoab
353: Collective
355: Input Parameters:
356: + dm - The DMMoab object being set
357: - range - The entities treated by this DMMoab
359: Level: beginner
361: @*/
362: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
363: {
364: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
367: dmmoab->elocal->clear();
368: dmmoab->eghost->clear();
369: dmmoab->elocal->insert(range->begin(), range->end());
370: #ifdef MOAB_HAVE_MPI
371: moab::ErrorCode merr;
372: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
373: MBERRNM(merr);
374: *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
375: #endif
376: dmmoab->neleloc = dmmoab->elocal->size();
377: dmmoab->neleghost = dmmoab->eghost->size();
378: #ifdef MOAB_HAVE_MPI
379: MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
380: PetscInfo(dm, "Created %" PetscInt_FMT " local and %" PetscInt_FMT " global elements.\n", dmmoab->neleloc, dmmoab->nele);
381: #else
382: dmmoab->nele = dmmoab->neleloc;
383: #endif
384: return 0;
385: }
387: /*@C
388: DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
390: Collective
392: Input Parameters:
393: + dm - The DMMoab object being set
394: - ltogtag - The MOAB tag used for local to global ids
396: Level: beginner
398: @*/
399: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
400: {
402: ((DM_Moab *)dm->data)->ltog_tag = ltogtag;
403: return 0;
404: }
406: /*@C
407: DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
409: Collective
411: Input Parameter:
412: . dm - The DMMoab object being set
414: Output Parameter:
415: . ltogtag - The MOAB tag used for local to global ids
417: Level: beginner
419: @*/
420: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
421: {
423: *ltog_tag = ((DM_Moab *)dm->data)->ltog_tag;
424: return 0;
425: }
427: /*@C
428: DMMoabSetBlockSize - Set the block size used with this DMMoab
430: Collective
432: Input Parameters:
433: + dm - The DMMoab object being set
434: - bs - The block size used with this DMMoab
436: Level: beginner
438: @*/
439: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
440: {
442: ((DM_Moab *)dm->data)->bs = bs;
443: return 0;
444: }
446: /*@C
447: DMMoabGetBlockSize - Get the block size used with this DMMoab
449: Collective
451: Input Parameter:
452: . dm - The DMMoab object being set
454: Output Parameter:
455: . bs - The block size used with this DMMoab
457: Level: beginner
459: @*/
460: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
461: {
463: *bs = ((DM_Moab *)dm->data)->bs;
464: return 0;
465: }
467: /*@C
468: DMMoabGetSize - Get the global vertex size used with this DMMoab
470: Collective on dm
472: Input Parameter:
473: . dm - The DMMoab object being set
475: Output Parameters:
476: + neg - The number of global elements in the DMMoab instance
477: - nvg - The number of global vertices in the DMMoab instance
479: Level: beginner
481: @*/
482: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
483: {
485: if (neg) *neg = ((DM_Moab *)dm->data)->nele;
486: if (nvg) *nvg = ((DM_Moab *)dm->data)->n;
487: return 0;
488: }
490: /*@C
491: DMMoabGetLocalSize - Get the local and ghosted vertex size used with this DMMoab
493: Collective on dm
495: Input Parameter:
496: . dm - The DMMoab object being set
498: Output Parameters:
499: + nel - The number of owned elements in this processor
500: . neg - The number of ghosted elements in this processor
501: . nvl - The number of owned vertices in this processor
502: - nvg - The number of ghosted vertices in this processor
504: Level: beginner
506: @*/
507: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
508: {
510: if (nel) *nel = ((DM_Moab *)dm->data)->neleloc;
511: if (neg) *neg = ((DM_Moab *)dm->data)->neleghost;
512: if (nvl) *nvl = ((DM_Moab *)dm->data)->nloc;
513: if (nvg) *nvg = ((DM_Moab *)dm->data)->nghost;
514: return 0;
515: }
517: /*@C
518: DMMoabGetOffset - Get the local offset for the global vector
520: Collective
522: Input Parameter:
523: . dm - The DMMoab object being set
525: Output Parameter:
526: . offset - The local offset for the global vector
528: Level: beginner
530: @*/
531: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
532: {
534: *offset = ((DM_Moab *)dm->data)->vstart;
535: return 0;
536: }
538: /*@C
539: DMMoabGetDimension - Get the dimension of the DM Mesh
541: Collective
543: Input Parameter:
544: . dm - The DMMoab object
546: Output Parameter:
547: . dim - The dimension of DM
549: Level: beginner
551: @*/
552: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
553: {
555: *dim = ((DM_Moab *)dm->data)->dim;
556: return 0;
557: }
559: /*@C
560: DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
561: generated through uniform refinement.
563: Collective on dm
565: Input Parameter:
566: . dm - The DMMoab object being set
568: Output Parameter:
569: . nvg - The current mesh hierarchy level
571: Level: beginner
573: @*/
574: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
575: {
577: if (nlevel) *nlevel = ((DM_Moab *)dm->data)->hlevel;
578: return 0;
579: }
581: /*@C
582: DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh
584: Collective
586: Input Parameters:
587: + dm - The DMMoab object
588: - ehandle - The element entity handle
590: Output Parameter:
591: . mat - The material ID for the current entity
593: Level: beginner
595: @*/
596: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
597: {
598: DM_Moab *dmmoab;
601: if (*mat) {
602: dmmoab = (DM_Moab *)(dm)->data;
603: *mat = dmmoab->materials[dmmoab->elocal->index(ehandle)];
604: }
605: return 0;
606: }
608: /*@C
609: DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities
611: Collective
613: Input Parameters:
614: + dm - The DMMoab object
615: . nconn - Number of entities whose coordinates are needed
616: - conn - The vertex entity handles
618: Output Parameter:
619: . vpos - The coordinates of the requested vertex entities
621: Level: beginner
623: .seealso: `DMMoabGetVertexConnectivity()`
624: @*/
625: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
626: {
627: DM_Moab *dmmoab;
628: moab::ErrorCode merr;
633: dmmoab = (DM_Moab *)(dm)->data;
635: /* Get connectivity information in MOAB canonical ordering */
636: if (dmmoab->hlevel) {
637: merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle *>(conn), nconn, dmmoab->hlevel, vpos);
638: MBERRNM(merr);
639: } else {
640: merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);
641: MBERRNM(merr);
642: }
643: return 0;
644: }
646: /*@C
647: DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity
649: Collective
651: Input Parameters:
652: + dm - The DMMoab object
653: - vhandle - Vertex entity handle
655: Output Parameters:
656: + nconn - Number of entities whose coordinates are needed
657: - conn - The vertex entity handles
659: Level: beginner
661: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabRestoreVertexConnectivity()`
662: @*/
663: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt *nconn, moab::EntityHandle **conn)
664: {
665: DM_Moab *dmmoab;
666: std::vector<moab::EntityHandle> adj_entities, connect;
667: moab::ErrorCode merr;
671: dmmoab = (DM_Moab *)(dm)->data;
673: /* Get connectivity information in MOAB canonical ordering */
674: merr = dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION);
675: MBERRNM(merr);
676: merr = dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect);
677: MBERRNM(merr);
679: if (conn) {
680: PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn);
681: PetscArraycpy(*conn, &connect[0], connect.size());
682: }
683: if (nconn) *nconn = connect.size();
684: return 0;
685: }
687: /*@C
688: DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity
690: Collective
692: Input Parameters:
693: + dm - The DMMoab object
694: . vhandle - Vertex entity handle
695: . nconn - Number of entities whose coordinates are needed
696: - conn - The vertex entity handles
698: Level: beginner
700: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`
701: @*/
702: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, moab::EntityHandle **conn)
703: {
707: if (conn) PetscFree(*conn);
708: if (nconn) *nconn = 0;
709: return 0;
710: }
712: /*@C
713: DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity
715: Collective
717: Input Parameters:
718: + dm - The DMMoab object
719: - ehandle - Vertex entity handle
721: Output Parameters:
722: + nconn - Number of entities whose coordinates are needed
723: - conn - The vertex entity handles
725: Level: beginner
727: .seealso: `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`, `DMMoabRestoreVertexConnectivity()`
728: @*/
729: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, const moab::EntityHandle **conn)
730: {
731: DM_Moab *dmmoab;
732: const moab::EntityHandle *connect;
733: std::vector<moab::EntityHandle> vconn;
734: moab::ErrorCode merr;
735: PetscInt nnodes;
739: dmmoab = (DM_Moab *)(dm)->data;
741: /* Get connectivity information in MOAB canonical ordering */
742: merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes);
743: MBERRNM(merr);
744: if (conn) *conn = connect;
745: if (nconn) *nconn = nnodes;
746: return 0;
747: }
749: /*@C
750: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
752: Collective
754: Input Parameters:
755: + dm - The DMMoab object
756: - ent - Entity handle
758: Output Parameter:
759: . ent_on_boundary - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
761: Level: beginner
763: .seealso: `DMMoabCheckBoundaryVertices()`
764: @*/
765: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool *ent_on_boundary)
766: {
767: moab::EntityType etype;
768: DM_Moab *dmmoab;
769: PetscInt edim;
773: dmmoab = (DM_Moab *)(dm)->data;
775: /* get the entity type and handle accordingly */
776: etype = dmmoab->mbiface->type_from_handle(ent);
779: /* get the entity dimension */
780: edim = dmmoab->mbiface->dimension_from_handle(ent);
782: *ent_on_boundary = PETSC_FALSE;
783: if (etype == moab::MBVERTEX && edim == 0) {
784: *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
785: } else {
786: if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
787: if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
788: } else { /* next check the lower-dimensional faces */
789: if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
790: }
791: }
792: return 0;
793: }
795: /*@C
796: DMMoabCheckBoundaryVertices - Check whether a given entity is on the boundary (vertex, edge, face, element)
798: Input Parameters:
799: + dm - The DMMoab object
800: . nconn - Number of handles
801: - cnt - Array of entity handles
803: Output Parameter:
804: . isbdvtx - Array of boundary markers - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
806: Level: beginner
808: .seealso: `DMMoabIsEntityOnBoundary()`
809: @*/
810: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool *isbdvtx)
811: {
812: DM_Moab *dmmoab;
813: PetscInt i;
818: dmmoab = (DM_Moab *)(dm)->data;
820: for (i = 0; i < nconn; ++i) isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
821: return 0;
822: }
824: /*@C
825: DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary
827: Input Parameter:
828: . dm - The DMMoab object
830: Output Parameters:
831: + bdvtx - Boundary vertices
832: . bdelems - Boundary elements
833: - bdfaces - Boundary faces
835: Level: beginner
837: .seealso: `DMMoabCheckBoundaryVertices()`, `DMMoabIsEntityOnBoundary()`
838: @*/
839: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range **bdelems, const moab::Range **bdfaces)
840: {
841: DM_Moab *dmmoab;
844: dmmoab = (DM_Moab *)(dm)->data;
846: if (bdvtx) *bdvtx = dmmoab->bndyvtx;
847: if (bdfaces) *bdfaces = dmmoab->bndyfaces;
848: if (bdelems) *bdfaces = dmmoab->bndyelems;
849: return 0;
850: }
852: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
853: {
854: PetscInt i;
855: moab::ErrorCode merr;
856: DM_Moab *dmmoab = (DM_Moab *)dm->data;
860: dmmoab->refct--;
861: if (!dmmoab->refct) {
862: delete dmmoab->vlocal;
863: delete dmmoab->vowned;
864: delete dmmoab->vghost;
865: delete dmmoab->elocal;
866: delete dmmoab->eghost;
867: delete dmmoab->bndyvtx;
868: delete dmmoab->bndyfaces;
869: delete dmmoab->bndyelems;
871: PetscFree(dmmoab->gsindices);
872: PetscFree2(dmmoab->gidmap, dmmoab->lidmap);
873: PetscFree(dmmoab->dfill);
874: PetscFree(dmmoab->ofill);
875: PetscFree(dmmoab->materials);
876: if (dmmoab->fieldNames) {
877: for (i = 0; i < dmmoab->numFields; i++) PetscFree(dmmoab->fieldNames[i]);
878: PetscFree(dmmoab->fieldNames);
879: }
881: if (dmmoab->nhlevels) {
882: PetscFree(dmmoab->hsets);
883: dmmoab->nhlevels = 0;
884: if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
885: dmmoab->hierarchy = NULL;
886: }
888: if (dmmoab->icreatedinstance) {
889: delete dmmoab->pcomm;
890: merr = dmmoab->mbiface->delete_mesh();
891: MBERRNM(merr);
892: delete dmmoab->mbiface;
893: }
894: dmmoab->mbiface = NULL;
895: #ifdef MOAB_HAVE_MPI
896: dmmoab->pcomm = NULL;
897: #endif
898: VecScatterDestroy(&dmmoab->ltog_sendrecv);
899: ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map);
900: PetscFree(dm->data);
901: }
902: return 0;
903: }
905: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(DM dm, PetscOptionItems *PetscOptionsObject)
906: {
907: DM_Moab *dmmoab = (DM_Moab *)dm->data;
909: PetscOptionsHeadBegin(PetscOptionsObject, "DMMoab Options");
910: PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL, 0);
911: 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);
912: /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
913: 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);
914: 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);
915: PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum *)&dmmoab->read_mode, NULL);
916: PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum *)&dmmoab->write_mode, NULL);
917: PetscOptionsHeadEnd();
918: return 0;
919: }
921: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
922: {
923: moab::ErrorCode merr;
924: Vec local, global;
925: IS from, to;
926: moab::Range::iterator iter;
927: PetscInt i, j, f, bs, vent, totsize, *lgmap;
928: DM_Moab *dmmoab = (DM_Moab *)dm->data;
929: moab::Range adjs;
932: /* Get the local and shared vertices and cache it */
934: #ifdef MOAB_HAVE_MPI
936: #endif
938: /* Get the entities recursively in the current part of the mesh, if user did not set the local vertices explicitly */
939: if (dmmoab->vlocal->empty()) {
940: //merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,*dmmoab->vlocal,true);MBERRNM(merr);
941: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, 0, *dmmoab->vlocal, false);
942: MBERRNM(merr);
944: #ifdef MOAB_HAVE_MPI
945: /* filter based on parallel status */
946: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
947: MBERRNM(merr);
949: /* filter all the non-owned and shared entities out of the list */
950: // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
951: adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
952: merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
953: MBERRNM(merr);
954: adjs = moab::subtract(adjs, *dmmoab->vghost);
955: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
956: #else
957: *dmmoab->vowned = *dmmoab->vlocal;
958: #endif
960: /* compute and cache the sizes of local and ghosted entities */
961: dmmoab->nloc = dmmoab->vowned->size();
962: dmmoab->nghost = dmmoab->vghost->size();
964: #ifdef MOAB_HAVE_MPI
965: MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
966: PetscInfo(NULL, "Filset ID: %lu, Vertices: local - %zu, owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->fileset, dmmoab->vlocal->size(), dmmoab->nloc, dmmoab->nghost);
967: #else
968: dmmoab->n = dmmoab->nloc;
969: #endif
970: }
972: {
973: /* get the information about the local elements in the mesh */
974: dmmoab->eghost->clear();
976: /* first decipher the leading dimension */
977: for (i = 3; i > 0; i--) {
978: dmmoab->elocal->clear();
979: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false);
980: MBERRNM(merr);
982: /* store the current mesh dimension */
983: if (dmmoab->elocal->size()) {
984: dmmoab->dim = i;
985: break;
986: }
987: }
989: DMSetDimension(dm, dmmoab->dim);
991: #ifdef MOAB_HAVE_MPI
992: /* filter the ghosted and owned element list */
993: *dmmoab->eghost = *dmmoab->elocal;
994: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
995: MBERRNM(merr);
996: *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
997: #endif
999: dmmoab->neleloc = dmmoab->elocal->size();
1000: dmmoab->neleghost = dmmoab->eghost->size();
1002: #ifdef MOAB_HAVE_MPI
1003: MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
1004: PetscInfo(NULL, "%d-dim elements: owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost);
1005: #else
1006: dmmoab->nele = dmmoab->neleloc;
1007: #endif
1008: }
1010: bs = dmmoab->bs;
1011: if (!dmmoab->ltog_tag) {
1012: /* Get the global ID tag. The global ID tag is applied to each
1013: vertex. It acts as an global identifier which MOAB uses to
1014: assemble the individual pieces of the mesh */
1015: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
1016: MBERRNM(merr);
1017: }
1019: totsize = dmmoab->vlocal->size();
1021: PetscCalloc1(totsize, &dmmoab->gsindices);
1022: {
1023: /* first get the local indices */
1024: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]);
1025: MBERRNM(merr);
1026: if (dmmoab->nghost) { /* next get the ghosted indices */
1027: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]);
1028: MBERRNM(merr);
1029: }
1031: /* find out the local and global minima of GLOBAL_ID */
1032: dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1033: for (i = 0; i < totsize; ++i) {
1034: if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1035: if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1036: }
1038: MPIU_Allreduce(&dmmoab->lminmax[0], &dmmoab->gminmax[0], 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm);
1039: MPIU_Allreduce(&dmmoab->lminmax[1], &dmmoab->gminmax[1], 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm);
1041: /* set the GID map */
1042: for (i = 0; i < totsize; ++i) { dmmoab->gsindices[i] -= dmmoab->gminmax[0]; /* zero based index needed for IS */ }
1043: dmmoab->lminmax[0] -= dmmoab->gminmax[0];
1044: dmmoab->lminmax[1] -= dmmoab->gminmax[0];
1046: 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]);
1047: }
1049: dmmoab->numFields);
1051: {
1052: dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1053: dmmoab->seqend = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1054: PetscInfo(NULL, "SEQUENCE: Local [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "]\n", dmmoab->seqstart, dmmoab->seqend);
1056: PetscMalloc2(dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->gidmap, dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->lidmap);
1057: PetscMalloc1(totsize * dmmoab->numFields, &lgmap);
1059: i = j = 0;
1060: /* set the owned vertex data first */
1061: for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1062: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1063: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1064: dmmoab->lidmap[vent] = i;
1065: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1066: }
1067: /* next arrange all the ghosted data information */
1068: for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1069: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1070: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1071: dmmoab->lidmap[vent] = i;
1072: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1073: }
1075: /* We need to create the Global to Local Vector Scatter Contexts
1076: 1) First create a local and global vector
1077: 2) Create a local and global IS
1078: 3) Create VecScatter and LtoGMapping objects
1079: 4) Cleanup the IS and Vec objects
1080: */
1081: DMCreateGlobalVector(dm, &global);
1082: DMCreateLocalVector(dm, &local);
1084: VecGetOwnershipRange(global, &dmmoab->vstart, &dmmoab->vend);
1086: /* global to local must retrieve ghost points */
1087: ISCreateStride(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, dmmoab->vstart, 1, &from);
1088: ISSetBlockSize(from, bs);
1090: ISCreateGeneral(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, &lgmap[0], PETSC_COPY_VALUES, &to);
1091: ISSetBlockSize(to, bs);
1093: if (!dmmoab->ltog_map) {
1094: /* create to the local to global mapping for vectors in order to use VecSetValuesLocal */
1095: ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm, dmmoab->bs, totsize * dmmoab->numFields, lgmap, PETSC_COPY_VALUES, &dmmoab->ltog_map);
1096: }
1098: /* now create the scatter object from local to global vector */
1099: VecScatterCreate(local, from, global, to, &dmmoab->ltog_sendrecv);
1101: /* clean up IS, Vec */
1102: PetscFree(lgmap);
1103: ISDestroy(&from);
1104: ISDestroy(&to);
1105: VecDestroy(&local);
1106: VecDestroy(&global);
1107: }
1109: dmmoab->bndyvtx = new moab::Range();
1110: dmmoab->bndyfaces = new moab::Range();
1111: dmmoab->bndyelems = new moab::Range();
1112: /* skin the boundary and store nodes */
1113: if (!dmmoab->hlevel) {
1114: /* get the skin vertices of boundary faces for the current partition and then filter
1115: the local, boundary faces, vertices and elements alone via PSTATUS flags;
1116: this should not give us any ghosted boundary, but if user needs such a functionality
1117: it would be easy to add it based on the find_skin query below */
1118: moab::Skinner skinner(dmmoab->mbiface);
1120: /* get the entities on the skin - only the faces */
1121: merr = skinner.find_skin(dmmoab->fileset, *dmmoab->elocal, false, *dmmoab->bndyfaces, NULL, true, true, false);
1122: MBERRNM(merr); // 'false' param indicates we want faces back, not vertices
1124: #ifdef MOAB_HAVE_MPI
1125: /* filter all the non-owned and shared entities out of the list */
1126: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1127: MBERRNM(merr);
1128: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT);
1129: MBERRNM(merr);
1130: #endif
1132: /* get all the nodes via connectivity and the parent elements via adjacency information */
1133: merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false);
1134: MBERRNM(merr);
1135: merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION);
1136: MBERRNM(merr);
1137: } else {
1138: /* Let us query the hierarchy manager and get the results directly for this level */
1139: for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1140: moab::EntityHandle elemHandle = *iter;
1141: if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1142: dmmoab->bndyelems->insert(elemHandle);
1143: /* For this boundary element, query the vertices and add them to the list */
1144: std::vector<moab::EntityHandle> connect;
1145: merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect);
1146: MBERRNM(merr);
1147: for (unsigned iv = 0; iv < connect.size(); ++iv)
1148: if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv])) dmmoab->bndyvtx->insert(connect[iv]);
1149: /* Next, let us query the boundary faces and add them also to the list */
1150: std::vector<moab::EntityHandle> faces;
1151: merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim - 1, faces);
1152: MBERRNM(merr);
1153: for (unsigned ifa = 0; ifa < faces.size(); ++ifa)
1154: if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa])) dmmoab->bndyfaces->insert(faces[ifa]);
1155: }
1156: }
1157: #ifdef MOAB_HAVE_MPI
1158: /* filter all the non-owned and shared entities out of the list */
1159: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1160: MBERRNM(merr);
1161: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1162: MBERRNM(merr);
1163: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1164: MBERRNM(merr);
1165: #endif
1166: }
1167: PetscInfo(NULL, "Found %zu boundary vertices, %zu boundary faces and %zu boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size());
1169: /* Get the material sets and populate the data for all locally owned elements */
1170: {
1171: PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials);
1172: /* Get the count of entities of particular type from dmmoab->elocal
1173: -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1174: moab::Range msets;
1175: merr = dmmoab->mbiface->get_entities_by_type_and_tag(dmmoab->fileset, moab::MBENTITYSET, &dmmoab->material_tag, NULL, 1, msets, moab::Interface::UNION);
1176: MB_CHK_ERR(merr);
1177: if (msets.size() == 0) PetscInfo(NULL, "No material sets found in the fileset.");
1179: for (unsigned i = 0; i < msets.size(); ++i) {
1180: moab::Range msetelems;
1181: merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);
1182: MB_CHK_ERR(merr);
1183: #ifdef MOAB_HAVE_MPI
1184: /* filter all the non-owned and shared entities out of the list */
1185: merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1186: MBERRNM(merr);
1187: #endif
1189: int partID;
1190: moab::EntityHandle mset = msets[i];
1191: merr = dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID);
1192: MB_CHK_ERR(merr);
1194: for (unsigned j = 0; j < msetelems.size(); ++j) dmmoab->materials[dmmoab->elocal->index(msetelems[j])] = partID;
1195: }
1196: }
1198: return 0;
1199: }
1201: /*@C
1202: DMMoabCreateVertices - Creates and adds several vertices to the primary set represented by the DM.
1204: Collective
1206: Input Parameters:
1207: + dm - The DM object
1208: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1209: . conn - The connectivity of the element
1210: - nverts - The number of vertices that form the element
1212: Output Parameter:
1213: . overts - The list of vertices that were created (can be NULL)
1215: Level: beginner
1217: .seealso: `DMMoabCreateSubmesh()`, `DMMoabCreateElement()`
1218: @*/
1219: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal *coords, PetscInt nverts, moab::Range *overts)
1220: {
1221: moab::ErrorCode merr;
1222: DM_Moab *dmmoab;
1223: moab::Range verts;
1228: dmmoab = (DM_Moab *)dm->data;
1230: /* Insert new points */
1231: merr = dmmoab->mbiface->create_vertices(&coords[0], nverts, verts);
1232: MBERRNM(merr);
1233: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, verts);
1234: MBERRNM(merr);
1236: if (overts) *overts = verts;
1237: return 0;
1238: }
1240: /*@C
1241: DMMoabCreateElement - Adds an element of specified type to the primary set represented by the DM.
1243: Collective
1245: Input Parameters:
1246: + dm - The DM object
1247: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1248: . conn - The connectivity of the element
1249: - nverts - The number of vertices that form the element
1251: Output Parameter:
1252: . oelem - The handle to the element created and added to the DM object
1254: Level: beginner
1256: .seealso: `DMMoabCreateSubmesh()`, `DMMoabCreateVertices()`
1257: @*/
1258: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle *conn, PetscInt nverts, moab::EntityHandle *oelem)
1259: {
1260: moab::ErrorCode merr;
1261: DM_Moab *dmmoab;
1262: moab::EntityHandle elem;
1267: dmmoab = (DM_Moab *)dm->data;
1269: /* Insert new element */
1270: merr = dmmoab->mbiface->create_element(type, conn, nverts, elem);
1271: MBERRNM(merr);
1272: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1);
1273: MBERRNM(merr);
1275: if (oelem) *oelem = elem;
1276: return 0;
1277: }
1279: /*@C
1280: DMMoabCreateSubmesh - Creates a sub-DM object with a set that contains all vertices/elements of the parent
1281: in addition to providing support for dynamic mesh modifications. This is useful for AMR calculations to
1282: create a DM object on a refined level.
1284: Collective
1286: Input Parameters:
1287: . dm - The DM object
1289: Output Parameter:
1290: . newdm - The sub DM object with updated set information
1292: Level: advanced
1294: .seealso: `DMCreate()`, `DMMoabCreateVertices()`, `DMMoabCreateElement()`
1295: @*/
1296: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1297: {
1298: DM_Moab *dmmoab;
1299: DM_Moab *ndmmoab;
1300: moab::ErrorCode merr;
1304: dmmoab = (DM_Moab *)dm->data;
1306: /* Create the basic DMMoab object and keep the default parameters created by DM impls */
1307: DMMoabCreateMoab(((PetscObject)dm)->comm, dmmoab->mbiface, &dmmoab->ltog_tag, PETSC_NULL, newdm);
1309: /* get all the necessary handles from the private DM object */
1310: ndmmoab = (DM_Moab *)(*newdm)->data;
1312: /* set the sub-mesh's parent DM reference */
1313: ndmmoab->parent = &dm;
1315: /* create a file set to associate all entities in current mesh */
1316: merr = ndmmoab->mbiface->create_meshset(moab::MESHSET_SET, ndmmoab->fileset);
1317: MBERR("Creating file set failed", merr);
1319: /* create a meshset and then add old fileset as child */
1320: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal);
1321: MBERR("Adding child vertices to parent failed", merr);
1322: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal);
1323: MBERR("Adding child elements to parent failed", merr);
1325: /* preserve the field association between the parent and sub-mesh objects */
1326: DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames);
1327: return 0;
1328: }
1330: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1331: {
1332: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
1333: const char *name;
1334: MPI_Comm comm;
1335: PetscMPIInt size;
1337: PetscObjectGetComm((PetscObject)dm, &comm);
1338: MPI_Comm_size(comm, &size);
1339: PetscObjectGetName((PetscObject)dm, &name);
1340: PetscViewerASCIIPushTab(viewer);
1341: if (name) PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimensions:\n", name, dmmoab->dim);
1342: else PetscViewerASCIIPrintf(viewer, "Mesh in %" PetscInt_FMT " dimensions:\n", dmmoab->dim);
1343: /* print details about the global mesh */
1344: {
1345: PetscViewerASCIIPushTab(viewer);
1346: PetscViewerASCIIPrintf(viewer, "Sizes: cells=%" PetscInt_FMT ", vertices=%" PetscInt_FMT ", blocks=%" PetscInt_FMT "\n", dmmoab->nele, dmmoab->n, dmmoab->bs);
1347: /* print boundary data */
1348: PetscViewerASCIIPrintf(viewer, "Boundary trace:\n");
1349: {
1350: PetscViewerASCIIPushTab(viewer);
1351: PetscViewerASCIIPrintf(viewer, "cells=%zu, faces=%zu, vertices=%zu\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size());
1352: PetscViewerASCIIPopTab(viewer);
1353: }
1354: /* print field data */
1355: PetscViewerASCIIPrintf(viewer, "Fields: %" PetscInt_FMT " components\n", dmmoab->numFields);
1356: {
1357: PetscViewerASCIIPushTab(viewer);
1358: for (int i = 0; i < dmmoab->numFields; ++i) PetscViewerASCIIPrintf(viewer, "[%" PetscInt_FMT "] - %s\n", i, dmmoab->fieldNames[i]);
1359: PetscViewerASCIIPopTab(viewer);
1360: }
1361: PetscViewerASCIIPopTab(viewer);
1362: }
1363: PetscViewerASCIIPopTab(viewer);
1364: PetscViewerFlush(viewer);
1365: return 0;
1366: }
1368: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1369: {
1370: return 0;
1371: }
1373: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1374: {
1375: return 0;
1376: }
1378: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1379: {
1380: PetscBool iascii, ishdf5, isvtk;
1384: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1385: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);
1386: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
1387: if (iascii) {
1388: DMMoabView_Ascii(dm, viewer);
1389: } else if (ishdf5) {
1390: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1391: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
1392: DMMoabView_HDF5(dm, viewer);
1393: PetscViewerPopFormat(viewer);
1394: #else
1395: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1396: #endif
1397: } else if (isvtk) {
1398: DMMoabView_VTK(dm, viewer);
1399: }
1400: return 0;
1401: }
1403: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
1404: {
1405: dm->ops->view = DMView_Moab;
1406: dm->ops->load = NULL /* DMLoad_Moab */;
1407: dm->ops->setfromoptions = DMSetFromOptions_Moab;
1408: dm->ops->clone = DMClone_Moab;
1409: dm->ops->setup = DMSetUp_Moab;
1410: dm->ops->createlocalsection = NULL;
1411: dm->ops->createdefaultconstraints = NULL;
1412: dm->ops->createglobalvector = DMCreateGlobalVector_Moab;
1413: dm->ops->createlocalvector = DMCreateLocalVector_Moab;
1414: dm->ops->getlocaltoglobalmapping = NULL;
1415: dm->ops->createfieldis = NULL;
1416: dm->ops->createcoordinatedm = NULL /* DMCreateCoordinateDM_Moab */;
1417: dm->ops->getcoloring = NULL;
1418: dm->ops->creatematrix = DMCreateMatrix_Moab;
1419: dm->ops->createinterpolation = DMCreateInterpolation_Moab;
1420: dm->ops->createinjection = NULL /* DMCreateInjection_Moab */;
1421: dm->ops->refine = DMRefine_Moab;
1422: dm->ops->coarsen = DMCoarsen_Moab;
1423: dm->ops->refinehierarchy = DMRefineHierarchy_Moab;
1424: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Moab;
1425: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Moab;
1426: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Moab;
1427: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Moab;
1428: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Moab;
1429: dm->ops->destroy = DMDestroy_Moab;
1430: dm->ops->createsubdm = NULL /* DMCreateSubDM_Moab */;
1431: dm->ops->getdimpoints = NULL /* DMGetDimPoints_Moab */;
1432: dm->ops->locatepoints = NULL /* DMLocatePoints_Moab */;
1433: return 0;
1434: }
1436: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1437: {
1438: /* get all the necessary handles from the private DM object */
1439: (*newdm)->data = (DM_Moab *)dm->data;
1440: ((DM_Moab *)dm->data)->refct++;
1442: PetscObjectChangeTypeName((PetscObject)*newdm, DMMOAB);
1443: DMInitialize_Moab(*newdm);
1444: return 0;
1445: }
1447: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1448: {
1450: PetscNew((DM_Moab **)&dm->data);
1452: ((DM_Moab *)dm->data)->bs = 1;
1453: ((DM_Moab *)dm->data)->numFields = 1;
1454: ((DM_Moab *)dm->data)->n = 0;
1455: ((DM_Moab *)dm->data)->nloc = 0;
1456: ((DM_Moab *)dm->data)->nghost = 0;
1457: ((DM_Moab *)dm->data)->nele = 0;
1458: ((DM_Moab *)dm->data)->neleloc = 0;
1459: ((DM_Moab *)dm->data)->neleghost = 0;
1460: ((DM_Moab *)dm->data)->ltog_map = NULL;
1461: ((DM_Moab *)dm->data)->ltog_sendrecv = NULL;
1463: ((DM_Moab *)dm->data)->refct = 1;
1464: ((DM_Moab *)dm->data)->parent = NULL;
1465: ((DM_Moab *)dm->data)->vlocal = new moab::Range();
1466: ((DM_Moab *)dm->data)->vowned = new moab::Range();
1467: ((DM_Moab *)dm->data)->vghost = new moab::Range();
1468: ((DM_Moab *)dm->data)->elocal = new moab::Range();
1469: ((DM_Moab *)dm->data)->eghost = new moab::Range();
1471: DMInitialize_Moab(dm);
1472: return 0;
1473: }