Actual source code: dmmoab.cxx
petsc-3.13.6 2020-09-29
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);
43: /* Un-implemented routines */
44: /*
45: PETSC_EXTERN PetscErrorCode DMCreatelocalsection_Moab(DM dm);
46: PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dmCoarse, DM dmFine, Mat *mat);
47: PETSC_EXTERN PetscErrorCode DMLoad_Moab(DM dm, PetscViewer viewer);
48: PETSC_EXTERN PetscErrorCode DMGetDimPoints_Moab(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd);
49: PETSC_EXTERN PetscErrorCode DMCreateSubDM_Moab(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
50: PETSC_EXTERN PetscErrorCode DMLocatePoints_Moab(DM dm, Vec v, IS *cellIS);
51: */
53: /*@C
54: DMMoabCreate - Creates a DMMoab object, which encapsulates a moab instance
56: Collective
58: Input Parameter:
59: . comm - The communicator for the DMMoab object
61: Output Parameter:
62: . dmb - The DMMoab object
64: Level: beginner
66: @*/
67: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
68: {
73: DMCreate(comm, dmb);
74: DMSetType(*dmb, DMMOAB);
75: return(0);
76: }
78: /*@C
79: DMMoabCreateMoab - Creates a DMMoab object, optionally from an instance and other data
81: Collective
83: Input Parameter:
84: + comm - The communicator for the DMMoab object
85: . mbiface - (ptr to) the MOAB Instance; if passed in NULL, MOAB instance is created inside PETSc, and destroyed
86: along with the DMMoab
87: . pcomm - (ptr to) a ParallelComm; if NULL, creates one internally for the whole communicator
88: . ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
89: - range - If non-NULL, contains range of entities to which DOFs will be assigned
91: Output Parameter:
92: . dmb - The DMMoab object
94: Level: intermediate
96: @*/
97: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
98: {
100: moab::ErrorCode merr;
101: DM dmmb;
102: DM_Moab *dmmoab;
107: DMMoabCreate(comm, &dmmb);
108: dmmoab = (DM_Moab*)(dmmb)->data;
110: if (!mbiface) {
111: dmmoab->mbiface = new moab::Core();
112: dmmoab->icreatedinstance = PETSC_TRUE;
113: }
114: else {
115: dmmoab->mbiface = mbiface;
116: dmmoab->icreatedinstance = PETSC_FALSE;
117: }
119: /* by default the fileset = root set. This set stores the hierarchy of entities belonging to current DM */
120: dmmoab->fileset = 0;
121: dmmoab->hlevel = 0;
122: dmmoab->nghostrings = 0;
124: #ifdef MOAB_HAVE_MPI
125: moab::EntityHandle partnset;
127: /* Create root sets for each mesh. Then pass these
128: to the load_file functions to be populated. */
129: merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, partnset); MBERR("Creating partition set failed", merr);
131: /* Create the parallel communicator object with the partition handle associated with MOAB */
132: dmmoab->pcomm = moab::ParallelComm::get_pcomm(dmmoab->mbiface, partnset, &comm);
133: #endif
135: /* do the remaining initializations for DMMoab */
136: dmmoab->bs = 1;
137: dmmoab->numFields = 1;
138: PetscMalloc(dmmoab->numFields * sizeof(char*), &dmmoab->fieldNames);
139: PetscStrallocpy("DEFAULT", (char**) &dmmoab->fieldNames[0]);
140: dmmoab->rw_dbglevel = 0;
141: dmmoab->partition_by_rank = PETSC_FALSE;
142: dmmoab->extra_read_options[0] = '\0';
143: dmmoab->extra_write_options[0] = '\0';
144: dmmoab->read_mode = READ_PART;
145: dmmoab->write_mode = WRITE_PART;
147: /* set global ID tag handle */
148: if (ltog_tag && *ltog_tag) {
149: DMMoabSetLocalToGlobalTag(dmmb, *ltog_tag);
150: }
151: else {
152: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag); MBERRNM(merr);
153: if (ltog_tag) *ltog_tag = dmmoab->ltog_tag;
154: }
156: merr = dmmoab->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dmmoab->material_tag); MBERRNM(merr);
158: /* set the local range of entities (vertices) of interest */
159: if (range) {
160: DMMoabSetLocalVertices(dmmb, range);
161: }
162: *dmb = dmmb;
163: return(0);
164: }
167: #ifdef MOAB_HAVE_MPI
169: /*@C
170: DMMoabGetParallelComm - Get the ParallelComm used with this DMMoab
172: Collective
174: Input Parameter:
175: . dm - The DMMoab object being set
177: Output Parameter:
178: . pcomm - The ParallelComm for the DMMoab
180: Level: beginner
182: @*/
183: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
184: {
187: *pcomm = ((DM_Moab*)(dm)->data)->pcomm;
188: return(0);
189: }
191: #endif /* MOAB_HAVE_MPI */
194: /*@C
195: DMMoabSetInterface - Set the MOAB instance used with this DMMoab
197: Collective
199: Input Parameter:
200: + dm - The DMMoab object being set
201: - mbiface - The MOAB instance being set on this DMMoab
203: Level: beginner
205: @*/
206: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
207: {
208: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
213: #ifdef MOAB_HAVE_MPI
214: dmmoab->pcomm = NULL;
215: #endif
216: dmmoab->mbiface = mbiface;
217: dmmoab->icreatedinstance = PETSC_FALSE;
218: return(0);
219: }
222: /*@C
223: DMMoabGetInterface - Get the MOAB instance used with this DMMoab
225: Collective
227: Input Parameter:
228: . dm - The DMMoab object being set
230: Output Parameter:
231: . mbiface - The MOAB instance set on this DMMoab
233: Level: beginner
235: @*/
236: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
237: {
238: PetscErrorCode ierr;
239: static PetscBool cite = PETSC_FALSE;
243: 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, K. and Stimpson, C. and Ernst, C.},\n year = {2004}, note = {Report}\n}\n", &cite);
244: *mbiface = ((DM_Moab*)dm->data)->mbiface;
245: return(0);
246: }
249: /*@C
250: DMMoabSetLocalVertices - Set the entities having DOFs on this DMMoab
252: Collective
254: Input Parameter:
255: + dm - The DMMoab object being set
256: - range - The entities treated by this DMMoab
258: Level: beginner
260: @*/
261: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
262: {
263: moab::Range tmpvtxs;
264: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
268: dmmoab->vlocal->clear();
269: dmmoab->vowned->clear();
271: dmmoab->vlocal->insert(range->begin(), range->end());
273: #ifdef MOAB_HAVE_MPI
274: moab::ErrorCode merr;
275: /* filter based on parallel status */
276: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned); MBERRNM(merr);
278: /* filter all the non-owned and shared entities out of the list */
279: tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
280: merr = dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost); MBERRNM(merr);
281: tmpvtxs = moab::subtract(tmpvtxs, *dmmoab->vghost);
282: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
283: #else
284: *dmmoab->vowned = *dmmoab->vlocal;
285: #endif
287: /* compute and cache the sizes of local and ghosted entities */
288: dmmoab->nloc = dmmoab->vowned->size();
289: dmmoab->nghost = dmmoab->vghost->size();
290: #ifdef MOAB_HAVE_MPI
291: PetscErrorCode ierr;
292: MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
293: #else
294: dmmoab->n = dmmoab->nloc;
295: #endif
296: return(0);
297: }
300: /*@C
301: DMMoabGetAllVertices - Get the entities having DOFs on this DMMoab
303: Collective
305: Input Parameter:
306: . dm - The DMMoab object being set
308: Output Parameter:
309: . owned - The local vertex entities in this DMMoab = (owned+ghosted)
311: Level: beginner
313: @*/
314: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
315: {
318: if (local) *local = *((DM_Moab*)dm->data)->vlocal;
319: return(0);
320: }
324: /*@C
325: DMMoabGetLocalVertices - Get the entities having DOFs on this DMMoab
327: Collective
329: Input Parameter:
330: . dm - The DMMoab object being set
332: Output Parameters:
333: + owned - The owned vertex entities in this DMMoab
334: - ghost - The ghosted entities (non-owned) stored locally in this partition
336: Level: beginner
338: @*/
339: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
340: {
343: if (owned) *owned = ((DM_Moab*)dm->data)->vowned;
344: if (ghost) *ghost = ((DM_Moab*)dm->data)->vghost;
345: return(0);
346: }
348: /*@C
349: DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned
351: Collective
353: Input Parameter:
354: . dm - The DMMoab object being set
356: Output Parameter:
357: . range - The entities owned locally
359: Level: beginner
361: @*/
362: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
363: {
366: if (range) *range = ((DM_Moab*)dm->data)->elocal;
367: return(0);
368: }
371: /*@C
372: DMMoabSetLocalElements - Set the entities having DOFs on this DMMoab
374: Collective
376: Input Parameters:
377: + dm - The DMMoab object being set
378: - range - The entities treated by this DMMoab
380: Level: beginner
382: @*/
383: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
384: {
385: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
389: dmmoab->elocal->clear();
390: dmmoab->eghost->clear();
391: dmmoab->elocal->insert(range->begin(), range->end());
392: #ifdef MOAB_HAVE_MPI
393: moab::ErrorCode merr;
394: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
395: *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
396: #endif
397: dmmoab->neleloc = dmmoab->elocal->size();
398: dmmoab->neleghost = dmmoab->eghost->size();
399: #ifdef MOAB_HAVE_MPI
400: PetscErrorCode ierr;
401: MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
402: PetscInfo2(dm, "Created %D local and %D global elements.\n", dmmoab->neleloc, dmmoab->nele);
403: #else
404: dmmoab->nele = dmmoab->neleloc;
405: #endif
406: return(0);
407: }
410: /*@C
411: DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
413: Collective
415: Input Parameters:
416: + dm - The DMMoab object being set
417: - ltogtag - The MOAB tag used for local to global ids
419: Level: beginner
421: @*/
422: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
423: {
426: ((DM_Moab*)dm->data)->ltog_tag = ltogtag;
427: return(0);
428: }
431: /*@C
432: DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
434: Collective
436: Input Parameter:
437: . dm - The DMMoab object being set
439: Output Parameter:
440: . ltogtag - The MOAB tag used for local to global ids
442: Level: beginner
444: @*/
445: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
446: {
449: *ltog_tag = ((DM_Moab*)dm->data)->ltog_tag;
450: return(0);
451: }
454: /*@C
455: DMMoabSetBlockSize - Set the block size used with this DMMoab
457: Collective
459: Input Parameter:
460: + dm - The DMMoab object being set
461: - bs - The block size used with this DMMoab
463: Level: beginner
465: @*/
466: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
467: {
470: ((DM_Moab*)dm->data)->bs = bs;
471: return(0);
472: }
475: /*@C
476: DMMoabGetBlockSize - Get the block size used with this DMMoab
478: Collective
480: Input Parameter:
481: . dm - The DMMoab object being set
483: Output Parameter:
484: . bs - The block size used with this DMMoab
486: Level: beginner
488: @*/
489: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
490: {
493: *bs = ((DM_Moab*)dm->data)->bs;
494: return(0);
495: }
498: /*@C
499: DMMoabGetSize - Get the global vertex size used with this DMMoab
501: Collective on dm
503: Input Parameter:
504: . dm - The DMMoab object being set
506: Output Parameter:
507: + neg - The number of global elements in the DMMoab instance
508: - nvg - The number of global vertices in the DMMoab instance
510: Level: beginner
512: @*/
513: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
514: {
517: if (neg) *neg = ((DM_Moab*)dm->data)->nele;
518: if (nvg) *nvg = ((DM_Moab*)dm->data)->n;
519: return(0);
520: }
523: /*@C
524: DMMoabGetLocalSize - Get the local and ghosted vertex size used with this DMMoab
526: Collective on dm
528: Input Parameter:
529: . dm - The DMMoab object being set
531: Output Parameter:
532: + nel - The number of owned elements in this processor
533: . neg - The number of ghosted elements in this processor
534: . nvl - The number of owned vertices in this processor
535: - nvg - The number of ghosted vertices in this processor
537: Level: beginner
539: @*/
540: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
541: {
544: if (nel) *nel = ((DM_Moab*)dm->data)->neleloc;
545: if (neg) *neg = ((DM_Moab*)dm->data)->neleghost;
546: if (nvl) *nvl = ((DM_Moab*)dm->data)->nloc;
547: if (nvg) *nvg = ((DM_Moab*)dm->data)->nghost;
548: return(0);
549: }
552: /*@C
553: DMMoabGetOffset - Get the local offset for the global vector
555: Collective
557: Input Parameter:
558: . dm - The DMMoab object being set
560: Output Parameter:
561: . offset - The local offset for the global vector
563: Level: beginner
565: @*/
566: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
567: {
570: *offset = ((DM_Moab*)dm->data)->vstart;
571: return(0);
572: }
575: /*@C
576: DMMoabGetDimension - Get the dimension of the DM Mesh
578: Collective
580: Input Parameter:
581: . dm - The DMMoab object
583: Output Parameter:
584: . dim - The dimension of DM
586: Level: beginner
588: @*/
589: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
590: {
593: *dim = ((DM_Moab*)dm->data)->dim;
594: return(0);
595: }
598: /*@C
599: DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
600: generated through uniform refinement.
602: Collective on dm
604: Input Parameter:
605: . dm - The DMMoab object being set
607: Output Parameter:
608: . nvg - The current mesh hierarchy level
610: Level: beginner
612: @*/
613: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
614: {
617: if (nlevel) *nlevel = ((DM_Moab*)dm->data)->hlevel;
618: return(0);
619: }
622: /*@C
623: DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh
625: Collective
627: Input Parameter:
628: + dm - The DMMoab object
629: - ehandle - The element entity handle
631: Output Parameter:
632: . mat - The material ID for the current entity
634: Level: beginner
636: @*/
637: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
638: {
639: DM_Moab *dmmoab;
643: if (*mat) {
644: dmmoab = (DM_Moab*)(dm)->data;
645: *mat = dmmoab->materials[dmmoab->elocal->index(ehandle)];
646: }
647: return(0);
648: }
651: /*@C
652: DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities
654: Collective
656: Input Parameter:
657: + dm - The DMMoab object
658: . nconn - Number of entities whose coordinates are needed
659: - conn - The vertex entity handles
661: Output Parameter:
662: . vpos - The coordinates of the requested vertex entities
664: Level: beginner
666: .seealso: DMMoabGetVertexConnectivity()
667: @*/
668: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
669: {
670: DM_Moab *dmmoab;
671: moab::ErrorCode merr;
677: dmmoab = (DM_Moab*)(dm)->data;
679: /* Get connectivity information in MOAB canonical ordering */
680: if (dmmoab->hlevel) {
681: merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle*>(conn), nconn, dmmoab->hlevel, vpos);MBERRNM(merr);
682: }
683: else {
684: merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);MBERRNM(merr);
685: }
686: return(0);
687: }
690: /*@C
691: DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity
693: Collective
695: Input Parameter:
696: + dm - The DMMoab object
697: - vhandle - Vertex entity handle
699: Output Parameter:
700: + nconn - Number of entities whose coordinates are needed
701: - conn - The vertex entity handles
703: Level: beginner
705: .seealso: DMMoabGetVertexCoordinates(), DMMoabRestoreVertexConnectivity()
706: @*/
707: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt* nconn, moab::EntityHandle **conn)
708: {
709: DM_Moab *dmmoab;
710: std::vector<moab::EntityHandle> adj_entities, connect;
711: PetscErrorCode ierr;
712: moab::ErrorCode merr;
717: dmmoab = (DM_Moab*)(dm)->data;
719: /* Get connectivity information in MOAB canonical ordering */
720: merr = dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION); MBERRNM(merr);
721: merr = dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect); MBERRNM(merr);
723: if (conn) {
724: PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn);
725: PetscArraycpy(*conn, &connect[0], connect.size());
726: }
727: if (nconn) *nconn = connect.size();
728: return(0);
729: }
732: /*@C
733: DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity
735: Collective
737: Input Parameter:
738: + dm - The DMMoab object
739: . vhandle - Vertex entity handle
740: . nconn - Number of entities whose coordinates are needed
741: - conn - The vertex entity handles
743: Level: beginner
745: .seealso: DMMoabGetVertexCoordinates(), DMMoabGetVertexConnectivity()
746: @*/
747: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt* nconn, moab::EntityHandle **conn)
748: {
749: PetscErrorCode ierr;
755: if (conn) {
756: PetscFree(*conn);
757: }
758: if (nconn) *nconn = 0;
759: return(0);
760: }
763: /*@C
764: DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity
766: Collective
768: Input Parameter:
769: + dm - The DMMoab object
770: - ehandle - Vertex entity handle
772: Output Parameter:
773: + nconn - Number of entities whose coordinates are needed
774: - conn - The vertex entity handles
776: Level: beginner
778: .seealso: DMMoabGetVertexCoordinates(), DMMoabGetVertexConnectivity(), DMMoabRestoreVertexConnectivity()
779: @*/
780: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt* nconn, const moab::EntityHandle **conn)
781: {
782: DM_Moab *dmmoab;
783: const moab::EntityHandle *connect;
784: std::vector<moab::EntityHandle> vconn;
785: moab::ErrorCode merr;
786: PetscInt nnodes;
791: dmmoab = (DM_Moab*)(dm)->data;
793: /* Get connectivity information in MOAB canonical ordering */
794: merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes); MBERRNM(merr);
795: if (conn) *conn = connect;
796: if (nconn) *nconn = nnodes;
797: return(0);
798: }
801: /*@C
802: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
804: Collective
806: Input Parameter:
807: + dm - The DMMoab object
808: - ent - Entity handle
810: Output Parameter:
811: . ent_on_boundary - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
813: Level: beginner
815: .seealso: DMMoabCheckBoundaryVertices()
816: @*/
817: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool* ent_on_boundary)
818: {
819: moab::EntityType etype;
820: DM_Moab *dmmoab;
821: PetscInt edim;
826: dmmoab = (DM_Moab*)(dm)->data;
828: /* get the entity type and handle accordingly */
829: etype = dmmoab->mbiface->type_from_handle(ent);
830: if (etype >= moab::MBPOLYHEDRON) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Entity type on the boundary skin is invalid. EntityType = %D\n", etype);
832: /* get the entity dimension */
833: edim = dmmoab->mbiface->dimension_from_handle(ent);
835: *ent_on_boundary = PETSC_FALSE;
836: if (etype == moab::MBVERTEX && edim == 0) {
837: *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
838: }
839: else {
840: if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
841: if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
842: }
843: else { /* next check the lower-dimensional faces */
844: if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
845: }
846: }
847: return(0);
848: }
851: /*@C
852: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
854: Input Parameter:
855: + dm - The DMMoab object
856: . nconn - Number of handles
857: - cnt - Array of entity handles
859: Output Parameter:
860: . isbdvtx - Array of boundary markers - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
862: Level: beginner
864: .seealso: DMMoabIsEntityOnBoundary()
865: @*/
866: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool* isbdvtx)
867: {
868: DM_Moab *dmmoab;
869: PetscInt i;
875: dmmoab = (DM_Moab*)(dm)->data;
877: for (i = 0; i < nconn; ++i) {
878: isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
879: }
880: return(0);
881: }
884: /*@C
885: DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary
887: Input Parameter:
888: . dm - The DMMoab object
890: Output Parameter:
891: + bdvtx - Boundary vertices
892: . bdelems - Boundary elements
893: - bdfaces - Boundary faces
895: Level: beginner
897: .seealso: DMMoabCheckBoundaryVertices(), DMMoabIsEntityOnBoundary()
898: @*/
899: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range** bdelems, const moab::Range** bdfaces)
900: {
901: DM_Moab *dmmoab;
905: dmmoab = (DM_Moab*)(dm)->data;
907: if (bdvtx) *bdvtx = dmmoab->bndyvtx;
908: if (bdfaces) *bdfaces = dmmoab->bndyfaces;
909: if (bdelems) *bdfaces = dmmoab->bndyelems;
910: return(0);
911: }
914: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
915: {
916: PetscErrorCode ierr;
917: PetscInt i;
918: moab::ErrorCode merr;
919: DM_Moab *dmmoab = (DM_Moab*)dm->data;
924: dmmoab->refct--;
925: if (!dmmoab->refct) {
926: delete dmmoab->vlocal;
927: delete dmmoab->vowned;
928: delete dmmoab->vghost;
929: delete dmmoab->elocal;
930: delete dmmoab->eghost;
931: delete dmmoab->bndyvtx;
932: delete dmmoab->bndyfaces;
933: delete dmmoab->bndyelems;
935: PetscFree(dmmoab->gsindices);
936: PetscFree2(dmmoab->gidmap, dmmoab->lidmap);
937: PetscFree(dmmoab->dfill);
938: PetscFree(dmmoab->ofill);
939: PetscFree(dmmoab->materials);
940: if (dmmoab->fieldNames) {
941: for (i = 0; i < dmmoab->numFields; i++) {
942: PetscFree(dmmoab->fieldNames[i]);
943: }
944: PetscFree(dmmoab->fieldNames);
945: }
947: if (dmmoab->nhlevels) {
948: PetscFree(dmmoab->hsets);
949: dmmoab->nhlevels = 0;
950: if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
951: dmmoab->hierarchy = NULL;
952: }
954: if (dmmoab->icreatedinstance) {
955: delete dmmoab->pcomm;
956: merr = dmmoab->mbiface->delete_mesh(); MBERRNM(merr);
957: delete dmmoab->mbiface;
958: }
959: dmmoab->mbiface = NULL;
960: #ifdef MOAB_HAVE_MPI
961: dmmoab->pcomm = NULL;
962: #endif
963: VecScatterDestroy(&dmmoab->ltog_sendrecv);
964: ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map);
965: PetscFree(dm->data);
966: }
967: return(0);
968: }
971: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(PetscOptionItems *PetscOptionsObject, DM dm)
972: {
974: DM_Moab *dmmoab = (DM_Moab*)dm->data;
978: PetscOptionsHead(PetscOptionsObject, "DMMoab Options");
979: PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL,0);
980: 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);
981: /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
982: 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, PETSC_MAX_PATH_LEN, NULL);
983: 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, PETSC_MAX_PATH_LEN, NULL);
984: PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum*)&dmmoab->read_mode, NULL);
985: PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum*)&dmmoab->write_mode, NULL);
986: return(0);
987: }
990: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
991: {
992: PetscErrorCode ierr;
993: moab::ErrorCode merr;
994: Vec local, global;
995: IS from, to;
996: moab::Range::iterator iter;
997: PetscInt i, j, f, bs, vent, totsize, *lgmap;
998: DM_Moab *dmmoab = (DM_Moab*)dm->data;
999: moab::Range adjs;
1003: /* Get the local and shared vertices and cache it */
1004: if (dmmoab->mbiface == NULL) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB Interface before calling SetUp.");
1005: #ifdef MOAB_HAVE_MPI
1006: if (dmmoab->pcomm == NULL) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB ParallelComm object before calling SetUp.");
1007: #endif
1009: /* Get the entities recursively in the current part of the mesh, if user did not set the local vertices explicitly */
1010: if (dmmoab->vlocal->empty())
1011: {
1012: //merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,*dmmoab->vlocal,true);MBERRNM(merr);
1013: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, 0, *dmmoab->vlocal, false); MBERRNM(merr);
1015: #ifdef MOAB_HAVE_MPI
1016: /* filter based on parallel status */
1017: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned); MBERRNM(merr);
1019: /* filter all the non-owned and shared entities out of the list */
1020: // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
1021: adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
1022: merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost); MBERRNM(merr);
1023: adjs = moab::subtract(adjs, *dmmoab->vghost);
1024: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
1025: #else
1026: *dmmoab->vowned = *dmmoab->vlocal;
1027: #endif
1029: /* compute and cache the sizes of local and ghosted entities */
1030: dmmoab->nloc = dmmoab->vowned->size();
1031: dmmoab->nghost = dmmoab->vghost->size();
1033: #ifdef MOAB_HAVE_MPI
1034: MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
1035: PetscInfo4(NULL, "Filset ID: %u, Vertices: local - %D, owned - %D, ghosted - %D.\n", dmmoab->fileset, dmmoab->vlocal->size(), dmmoab->nloc, dmmoab->nghost);
1036: #else
1037: dmmoab->n = dmmoab->nloc;
1038: #endif
1039: }
1041: {
1042: /* get the information about the local elements in the mesh */
1043: dmmoab->eghost->clear();
1045: /* first decipher the leading dimension */
1046: for (i = 3; i > 0; i--) {
1047: dmmoab->elocal->clear();
1048: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false); MBERRNM(merr);
1050: /* store the current mesh dimension */
1051: if (dmmoab->elocal->size()) {
1052: dmmoab->dim = i;
1053: break;
1054: }
1055: }
1057: DMSetDimension(dm, dmmoab->dim);
1059: #ifdef MOAB_HAVE_MPI
1060: /* filter the ghosted and owned element list */
1061: *dmmoab->eghost = *dmmoab->elocal;
1062: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1063: *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
1064: #endif
1066: dmmoab->neleloc = dmmoab->elocal->size();
1067: dmmoab->neleghost = dmmoab->eghost->size();
1069: #ifdef MOAB_HAVE_MPI
1070: MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
1071: PetscInfo3(NULL, "%d-dim elements: owned - %D, ghosted - %D.\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost);
1072: #else
1073: dmmoab->nele = dmmoab->neleloc;
1074: #endif
1075: }
1077: bs = dmmoab->bs;
1078: if (!dmmoab->ltog_tag) {
1079: /* Get the global ID tag. The global ID tag is applied to each
1080: vertex. It acts as an global identifier which MOAB uses to
1081: assemble the individual pieces of the mesh */
1082: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag); MBERRNM(merr);
1083: }
1085: totsize = dmmoab->vlocal->size();
1086: if (totsize != dmmoab->nloc + dmmoab->nghost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Mismatch between local and owned+ghost vertices. %D != %D.", totsize, dmmoab->nloc + dmmoab->nghost);
1087: PetscCalloc1(totsize, &dmmoab->gsindices);
1088: {
1089: /* first get the local indices */
1090: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]); MBERRNM(merr);
1091: if (dmmoab->nghost) { /* next get the ghosted indices */
1092: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]); MBERRNM(merr);
1093: }
1095: /* find out the local and global minima of GLOBAL_ID */
1096: dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1097: for (i = 0; i < totsize; ++i) {
1098: if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1099: if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1100: }
1102: MPIU_Allreduce(&dmmoab->lminmax[0], &dmmoab->gminmax[0], 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm);
1103: MPIU_Allreduce(&dmmoab->lminmax[1], &dmmoab->gminmax[1], 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm);
1105: /* set the GID map */
1106: for (i = 0; i < totsize; ++i) {
1107: dmmoab->gsindices[i] -= dmmoab->gminmax[0]; /* zero based index needed for IS */
1109: }
1110: dmmoab->lminmax[0] -= dmmoab->gminmax[0];
1111: dmmoab->lminmax[1] -= dmmoab->gminmax[0];
1113: PetscInfo4(NULL, "GLOBAL_ID: Local [min, max] - [%D, %D], Global [min, max] - [%D, %D]\n", dmmoab->lminmax[0], dmmoab->lminmax[1], dmmoab->gminmax[0], dmmoab->gminmax[1]);
1114: }
1115: if (!(dmmoab->bs == dmmoab->numFields || dmmoab->bs == 1)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mismatch between block size and number of component fields. %D != 1 OR %D != %D.", dmmoab->bs, dmmoab->bs, dmmoab->numFields);
1117: {
1118: dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1119: dmmoab->seqend = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1120: PetscInfo2(NULL, "SEQUENCE: Local [min, max] - [%D, %D]\n", dmmoab->seqstart, dmmoab->seqend);
1122: PetscMalloc2(dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->gidmap, dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->lidmap);
1123: PetscMalloc1(totsize * dmmoab->numFields, &lgmap);
1125: i = j = 0;
1126: /* set the owned vertex data first */
1127: for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1128: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1129: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1130: dmmoab->lidmap[vent] = i;
1131: for (f = 0; f < dmmoab->numFields; f++, j++) {
1132: lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1133: }
1134: }
1135: /* next arrange all the ghosted data information */
1136: for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1137: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1138: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1139: dmmoab->lidmap[vent] = i;
1140: for (f = 0; f < dmmoab->numFields; f++, j++) {
1141: lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1142: }
1143: }
1145: /* We need to create the Global to Local Vector Scatter Contexts
1146: 1) First create a local and global vector
1147: 2) Create a local and global IS
1148: 3) Create VecScatter and LtoGMapping objects
1149: 4) Cleanup the IS and Vec objects
1150: */
1151: DMCreateGlobalVector(dm, &global);
1152: DMCreateLocalVector(dm, &local);
1154: VecGetOwnershipRange(global, &dmmoab->vstart, &dmmoab->vend);
1156: /* global to local must retrieve ghost points */
1157: ISCreateStride(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, dmmoab->vstart, 1, &from);
1158: ISSetBlockSize(from, bs);
1160: ISCreateGeneral(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, &lgmap[0], PETSC_COPY_VALUES, &to);
1161: ISSetBlockSize(to, bs);
1163: if (!dmmoab->ltog_map) {
1164: /* create to the local to global mapping for vectors in order to use VecSetValuesLocal */
1165: ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm, dmmoab->bs, totsize * dmmoab->numFields, lgmap,
1166: PETSC_COPY_VALUES, &dmmoab->ltog_map);
1167: }
1169: /* now create the scatter object from local to global vector */
1170: VecScatterCreate(local, from, global, to, &dmmoab->ltog_sendrecv);
1172: /* clean up IS, Vec */
1173: PetscFree(lgmap);
1174: ISDestroy(&from);
1175: ISDestroy(&to);
1176: VecDestroy(&local);
1177: VecDestroy(&global);
1178: }
1180: dmmoab->bndyvtx = new moab::Range();
1181: dmmoab->bndyfaces = new moab::Range();
1182: dmmoab->bndyelems = new moab::Range();
1183: /* skin the boundary and store nodes */
1184: if (!dmmoab->hlevel) {
1185: /* get the skin vertices of boundary faces for the current partition and then filter
1186: the local, boundary faces, vertices and elements alone via PSTATUS flags;
1187: this should not give us any ghosted boundary, but if user needs such a functionality
1188: it would be easy to add it based on the find_skin query below */
1189: moab::Skinner skinner(dmmoab->mbiface);
1191: /* get the entities on the skin - only the faces */
1192: merr = skinner.find_skin(dmmoab->fileset, *dmmoab->elocal, false, *dmmoab->bndyfaces, NULL, true, true, false); MBERRNM(merr); // 'false' param indicates we want faces back, not vertices
1194: #ifdef MOAB_HAVE_MPI
1195: /* filter all the non-owned and shared entities out of the list */
1196: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1197: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT); MBERRNM(merr);
1198: #endif
1200: /* get all the nodes via connectivity and the parent elements via adjacency information */
1201: merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false); MBERRNM(ierr);
1202: merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION); MBERRNM(ierr);
1203: }
1204: else {
1205: /* Let us query the hierarchy manager and get the results directly for this level */
1206: for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1207: moab::EntityHandle elemHandle = *iter;
1208: if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1209: dmmoab->bndyelems->insert(elemHandle);
1210: /* For this boundary element, query the vertices and add them to the list */
1211: std::vector<moab::EntityHandle> connect;
1212: merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect); MBERRNM(ierr);
1213: for (unsigned iv=0; iv < connect.size(); ++iv)
1214: if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv]))
1215: dmmoab->bndyvtx->insert(connect[iv]);
1216: /* Next, let us query the boundary faces and add them also to the list */
1217: std::vector<moab::EntityHandle> faces;
1218: merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim-1, faces); MBERRNM(ierr);
1219: for (unsigned ifa=0; ifa < faces.size(); ++ifa)
1220: if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa]))
1221: dmmoab->bndyfaces->insert(faces[ifa]);
1222: }
1223: }
1224: #ifdef MOAB_HAVE_MPI
1225: /* filter all the non-owned and shared entities out of the list */
1226: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1227: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1228: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1229: #endif
1231: }
1232: PetscInfo3(NULL, "Found %D boundary vertices, %D boundary faces and %D boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size());
1234: /* Get the material sets and populate the data for all locally owned elements */
1235: {
1236: PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials);
1237: /* Get the count of entities of particular type from dmmoab->elocal
1238: -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1239: moab::Range msets;
1240: merr = dmmoab->mbiface->get_entities_by_type_and_tag(dmmoab->fileset, moab::MBENTITYSET, &dmmoab->material_tag, NULL, 1, msets, moab::Interface::UNION);MB_CHK_ERR(merr);
1241: if (msets.size() == 0) {
1242: PetscInfo(NULL, "No material sets found in the fileset.");
1243: }
1245: for (unsigned i=0; i < msets.size(); ++i) {
1246: moab::Range msetelems;
1247: merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);MB_CHK_ERR(merr);
1248: #ifdef MOAB_HAVE_MPI
1249: /* filter all the non-owned and shared entities out of the list */
1250: merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1251: #endif
1253: int partID;
1254: moab::EntityHandle mset=msets[i];
1255: merr = dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID);MB_CHK_ERR(merr);
1257: for (unsigned j=0; j < msetelems.size(); ++j)
1258: dmmoab->materials[dmmoab->elocal->index(msetelems[j])]=partID;
1259: }
1260: }
1262: return(0);
1263: }
1266: /*@C
1267: DMMoabCreateVertices - Creates and adds several vertices to the primary set represented by the DM.
1269: Collective
1271: Input Parameters:
1272: + dm - The DM object
1273: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1274: . conn - The connectivity of the element
1275: - nverts - The number of vertices that form the element
1277: Output Parameter:
1278: . overts - The list of vertices that were created (can be NULL)
1280: Level: beginner
1282: .seealso: DMMoabCreateSubmesh(), DMMoabCreateElement()
1283: @*/
1284: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal* coords, PetscInt nverts, moab::Range* overts)
1285: {
1286: moab::ErrorCode merr;
1287: DM_Moab *dmmoab;
1288: moab::Range verts;
1294: dmmoab = (DM_Moab*) dm->data;
1296: /* Insert new points */
1297: merr = dmmoab->mbiface->create_vertices(&coords[0], nverts, verts); MBERRNM(merr);
1298: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, verts); MBERRNM(merr);
1300: if (overts) *overts = verts;
1301: return(0);
1302: }
1305: /*@C
1306: DMMoabCreateElement - Adds an element of specified type to the primary set represented by the DM.
1308: Collective
1310: Input Parameters:
1311: + dm - The DM object
1312: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1313: . conn - The connectivity of the element
1314: - nverts - The number of vertices that form the element
1316: Output Parameter:
1317: . oelem - The handle to the element created and added to the DM object
1319: Level: beginner
1321: .seealso: DMMoabCreateSubmesh(), DMMoabCreateVertices()
1322: @*/
1323: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle* conn, PetscInt nverts, moab::EntityHandle* oelem)
1324: {
1325: moab::ErrorCode merr;
1326: DM_Moab *dmmoab;
1327: moab::EntityHandle elem;
1333: dmmoab = (DM_Moab*) dm->data;
1335: /* Insert new element */
1336: merr = dmmoab->mbiface->create_element(type, conn, nverts, elem); MBERRNM(merr);
1337: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1); MBERRNM(merr);
1339: if (oelem) *oelem = elem;
1340: return(0);
1341: }
1344: /*@C
1345: DMMoabCreateSubmesh - Creates a sub-DM object with a set that contains all vertices/elements of the parent
1346: in addition to providing support for dynamic mesh modifications. This is useful for AMR calculations to
1347: create a DM object on a refined level.
1349: Collective
1351: Input Parameters:
1352: . dm - The DM object
1354: Output Parameter:
1355: . newdm - The sub DM object with updated set information
1357: Level: advanced
1359: .seealso: DMCreate(), DMMoabCreateVertices(), DMMoabCreateElement()
1360: @*/
1361: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1362: {
1363: DM_Moab *dmmoab;
1364: DM_Moab *ndmmoab;
1365: moab::ErrorCode merr;
1366: PetscErrorCode ierr;
1371: dmmoab = (DM_Moab*) dm->data;
1373: /* Create the basic DMMoab object and keep the default parameters created by DM impls */
1374: DMMoabCreateMoab(((PetscObject)dm)->comm, dmmoab->mbiface, &dmmoab->ltog_tag, PETSC_NULL, newdm);
1376: /* get all the necessary handles from the private DM object */
1377: ndmmoab = (DM_Moab*) (*newdm)->data;
1379: /* set the sub-mesh's parent DM reference */
1380: ndmmoab->parent = &dm;
1382: /* create a file set to associate all entities in current mesh */
1383: merr = ndmmoab->mbiface->create_meshset(moab::MESHSET_SET, ndmmoab->fileset); MBERR("Creating file set failed", merr);
1385: /* create a meshset and then add old fileset as child */
1386: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal); MBERR("Adding child vertices to parent failed", merr);
1387: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal); MBERR("Adding child elements to parent failed", merr);
1389: /* preserve the field association between the parent and sub-mesh objects */
1390: DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames);
1391: return(0);
1392: }
1395: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1396: {
1397: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
1398: const char *name;
1399: MPI_Comm comm;
1400: PetscMPIInt size;
1401: PetscErrorCode ierr;
1404: PetscObjectGetComm((PetscObject)dm, &comm);
1405: MPI_Comm_size(comm, &size);
1406: PetscObjectGetName((PetscObject) dm, &name);
1407: PetscViewerASCIIPushTab(viewer);
1408: if (name) {PetscViewerASCIIPrintf(viewer, "%s in %D dimensions:\n", name, dmmoab->dim);}
1409: else {PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dmmoab->dim);}
1410: /* print details about the global mesh */
1411: {
1412: PetscViewerASCIIPushTab(viewer);
1413: PetscViewerASCIIPrintf(viewer, "Sizes: cells=%D, vertices=%D, blocks=%D\n", dmmoab->nele, dmmoab->n, dmmoab->bs);
1414: /* print boundary data */
1415: PetscViewerASCIIPrintf(viewer, "Boundary trace:\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size());
1416: {
1417: PetscViewerASCIIPushTab(viewer);
1418: PetscViewerASCIIPrintf(viewer, "cells=%D, faces=%D, vertices=%D\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size());
1419: PetscViewerASCIIPopTab(viewer);
1420: }
1421: /* print field data */
1422: PetscViewerASCIIPrintf(viewer, "Fields: %D components\n", dmmoab->numFields);
1423: {
1424: PetscViewerASCIIPushTab(viewer);
1425: for (int i = 0; i < dmmoab->numFields; ++i) {
1426: PetscViewerASCIIPrintf(viewer, "[%D] - %s\n", i, dmmoab->fieldNames[i]);
1427: }
1428: PetscViewerASCIIPopTab(viewer);
1429: }
1430: PetscViewerASCIIPopTab(viewer);
1431: }
1432: PetscViewerASCIIPopTab(viewer);
1433: PetscViewerFlush(viewer);
1434: return(0);
1435: }
1437: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1438: {
1439: return(0);
1440: }
1442: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1443: {
1444: return(0);
1445: }
1447: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1448: {
1449: PetscBool iascii, ishdf5, isvtk;
1455: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1456: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
1457: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
1458: if (iascii) {
1459: DMMoabView_Ascii(dm, viewer);
1460: } else if (ishdf5) {
1461: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1462: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
1463: DMMoabView_HDF5(dm, viewer);
1464: PetscViewerPopFormat(viewer);
1465: #else
1466: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1467: #endif
1468: }
1469: else if (isvtk) {
1470: DMMoabView_VTK(dm, viewer);
1471: }
1472: return(0);
1473: }
1476: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
1477: {
1479: dm->ops->view = DMView_Moab;
1480: dm->ops->load = NULL /* DMLoad_Moab */;
1481: dm->ops->setfromoptions = DMSetFromOptions_Moab;
1482: dm->ops->clone = DMClone_Moab;
1483: dm->ops->setup = DMSetUp_Moab;
1484: dm->ops->createlocalsection = NULL;
1485: dm->ops->createdefaultconstraints = NULL;
1486: dm->ops->createglobalvector = DMCreateGlobalVector_Moab;
1487: dm->ops->createlocalvector = DMCreateLocalVector_Moab;
1488: dm->ops->getlocaltoglobalmapping = NULL;
1489: dm->ops->createfieldis = NULL;
1490: dm->ops->createcoordinatedm = NULL /* DMCreateCoordinateDM_Moab */;
1491: dm->ops->getcoloring = NULL;
1492: dm->ops->creatematrix = DMCreateMatrix_Moab;
1493: dm->ops->createinterpolation = DMCreateInterpolation_Moab;
1494: dm->ops->createinjection = NULL /* DMCreateInjection_Moab */;
1495: dm->ops->refine = DMRefine_Moab;
1496: dm->ops->coarsen = DMCoarsen_Moab;
1497: dm->ops->refinehierarchy = DMRefineHierarchy_Moab;
1498: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Moab;
1499: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Moab;
1500: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Moab;
1501: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Moab;
1502: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Moab;
1503: dm->ops->destroy = DMDestroy_Moab;
1504: dm->ops->createsubdm = NULL /* DMCreateSubDM_Moab */;
1505: dm->ops->getdimpoints = NULL /* DMGetDimPoints_Moab */;
1506: dm->ops->locatepoints = NULL /* DMLocatePoints_Moab */;
1507: return(0);
1508: }
1511: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1512: {
1513: PetscErrorCode ierr;
1516: /* get all the necessary handles from the private DM object */
1517: (*newdm)->data = (DM_Moab*) dm->data;
1518: ((DM_Moab*)dm->data)->refct++;
1520: PetscObjectChangeTypeName((PetscObject) *newdm, DMMOAB);
1521: DMInitialize_Moab(*newdm);
1522: return(0);
1523: }
1526: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1527: {
1532: PetscNewLog(dm, (DM_Moab**)&dm->data);
1534: ((DM_Moab*)dm->data)->bs = 1;
1535: ((DM_Moab*)dm->data)->numFields = 1;
1536: ((DM_Moab*)dm->data)->n = 0;
1537: ((DM_Moab*)dm->data)->nloc = 0;
1538: ((DM_Moab*)dm->data)->nghost = 0;
1539: ((DM_Moab*)dm->data)->nele = 0;
1540: ((DM_Moab*)dm->data)->neleloc = 0;
1541: ((DM_Moab*)dm->data)->neleghost = 0;
1542: ((DM_Moab*)dm->data)->ltog_map = NULL;
1543: ((DM_Moab*)dm->data)->ltog_sendrecv = NULL;
1545: ((DM_Moab*)dm->data)->refct = 1;
1546: ((DM_Moab*)dm->data)->parent = NULL;
1547: ((DM_Moab*)dm->data)->vlocal = new moab::Range();
1548: ((DM_Moab*)dm->data)->vowned = new moab::Range();
1549: ((DM_Moab*)dm->data)->vghost = new moab::Range();
1550: ((DM_Moab*)dm->data)->elocal = new moab::Range();
1551: ((DM_Moab*)dm->data)->eghost = new moab::Range();
1553: DMInitialize_Moab(dm);
1554: return(0);
1555: }