Actual source code: dmmoab.cxx
petsc-3.11.4 2019-09-28
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: http://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 DMCreateDefaultSection_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: /*@
54: DMMoabCreate - Creates a DMMoab object, which encapsulates a moab instance
56: Collective on MPI_Comm
58: Input Parameter:
59: . comm - The communicator for the DMMoab object
61: Output Parameter:
62: . dmb - The DMMoab object
64: Level: beginner
66: .keywords: DMMoab, create
67: @*/
68: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
69: {
74: DMCreate(comm, dmb);
75: DMSetType(*dmb, DMMOAB);
76: return(0);
77: }
79: /*@
80: DMMoabCreateMoab - Creates a DMMoab object, optionally from an instance and other data
82: Collective on MPI_Comm
84: Input Parameter:
85: . comm - The communicator for the DMMoab object
86: . mbiface - (ptr to) the MOAB Instance; if passed in NULL, MOAB instance is created inside PETSc, and destroyed
87: along with the DMMoab
88: . pcomm - (ptr to) a ParallelComm; if NULL, creates one internally for the whole communicator
89: . ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
90: . range - If non-NULL, contains range of entities to which DOFs will be assigned
92: Output Parameter:
93: . dmb - The DMMoab object
95: Level: intermediate
97: .keywords: DMMoab, create
98: @*/
99: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
100: {
102: moab::ErrorCode merr;
103: DM dmmb;
104: DM_Moab *dmmoab;
109: DMMoabCreate(comm, &dmmb);
110: dmmoab = (DM_Moab*)(dmmb)->data;
112: if (!mbiface) {
113: dmmoab->mbiface = new moab::Core();
114: dmmoab->icreatedinstance = PETSC_TRUE;
115: }
116: else {
117: dmmoab->mbiface = mbiface;
118: dmmoab->icreatedinstance = PETSC_FALSE;
119: }
121: /* by default the fileset = root set. This set stores the hierarchy of entities belonging to current DM */
122: dmmoab->fileset = 0;
123: dmmoab->hlevel = 0;
124: dmmoab->nghostrings = 0;
126: #ifdef MOAB_HAVE_MPI
127: moab::EntityHandle partnset;
129: /* Create root sets for each mesh. Then pass these
130: to the load_file functions to be populated. */
131: merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, partnset); MBERR("Creating partition set failed", merr);
133: /* Create the parallel communicator object with the partition handle associated with MOAB */
134: dmmoab->pcomm = moab::ParallelComm::get_pcomm(dmmoab->mbiface, partnset, &comm);
135: #endif
137: /* do the remaining initializations for DMMoab */
138: dmmoab->bs = 1;
139: dmmoab->numFields = 1;
140: PetscMalloc(dmmoab->numFields * sizeof(char*), &dmmoab->fieldNames);
141: PetscStrallocpy("DEFAULT", (char**) &dmmoab->fieldNames[0]);
142: dmmoab->rw_dbglevel = 0;
143: dmmoab->partition_by_rank = PETSC_FALSE;
144: dmmoab->extra_read_options[0] = '\0';
145: dmmoab->extra_write_options[0] = '\0';
146: dmmoab->read_mode = READ_PART;
147: dmmoab->write_mode = WRITE_PART;
149: /* set global ID tag handle */
150: if (ltog_tag && *ltog_tag) {
151: DMMoabSetLocalToGlobalTag(dmmb, *ltog_tag);
152: }
153: else {
154: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag); MBERRNM(merr);
155: if (ltog_tag) *ltog_tag = dmmoab->ltog_tag;
156: }
158: merr = dmmoab->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dmmoab->material_tag); MBERRNM(merr);
160: /* set the local range of entities (vertices) of interest */
161: if (range) {
162: DMMoabSetLocalVertices(dmmb, range);
163: }
164: *dmb = dmmb;
165: return(0);
166: }
169: #ifdef MOAB_HAVE_MPI
171: /*@
172: DMMoabGetParallelComm - Get the ParallelComm used with this DMMoab
174: Collective on MPI_Comm
176: Input Parameter:
177: . dm - The DMMoab object being set
179: Output Parameter:
180: . pcomm - The ParallelComm for the DMMoab
182: Level: beginner
184: .keywords: DMMoab, create
185: @*/
186: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
187: {
190: *pcomm = ((DM_Moab*)(dm)->data)->pcomm;
191: return(0);
192: }
194: #endif /* MOAB_HAVE_MPI */
197: /*@
198: DMMoabSetInterface - Set the MOAB instance used with this DMMoab
200: Collective on MPI_Comm
202: Input Parameter:
203: . dm - The DMMoab object being set
204: . mbiface - The MOAB instance being set on this DMMoab
206: Level: beginner
208: .keywords: DMMoab, create
209: @*/
210: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
211: {
212: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
217: #ifdef MOAB_HAVE_MPI
218: dmmoab->pcomm = NULL;
219: #endif
220: dmmoab->mbiface = mbiface;
221: dmmoab->icreatedinstance = PETSC_FALSE;
222: return(0);
223: }
226: /*@
227: DMMoabGetInterface - Get the MOAB instance used with this DMMoab
229: Collective on MPI_Comm
231: Input Parameter:
232: . dm - The DMMoab object being set
234: Output Parameter:
235: . mbiface - The MOAB instance set on this DMMoab
237: Level: beginner
239: .keywords: DMMoab, create
240: @*/
241: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
242: {
243: PetscErrorCode ierr;
244: static PetscBool cite = PETSC_FALSE;
248: 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);
249: *mbiface = ((DM_Moab*)dm->data)->mbiface;
250: return(0);
251: }
254: /*@
255: DMMoabSetLocalVertices - Set the entities having DOFs on this DMMoab
257: Collective on MPI_Comm
259: Input Parameter:
260: . dm - The DMMoab object being set
261: . range - The entities treated by this DMMoab
263: Level: beginner
265: .keywords: DMMoab, create
266: @*/
267: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
268: {
269: moab::Range tmpvtxs;
270: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
274: dmmoab->vlocal->clear();
275: dmmoab->vowned->clear();
277: dmmoab->vlocal->insert(range->begin(), range->end());
279: #ifdef MOAB_HAVE_MPI
280: moab::ErrorCode merr;
281: /* filter based on parallel status */
282: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned); MBERRNM(merr);
284: /* filter all the non-owned and shared entities out of the list */
285: tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
286: merr = dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost); MBERRNM(merr);
287: tmpvtxs = moab::subtract(tmpvtxs, *dmmoab->vghost);
288: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
289: #else
290: *dmmoab->vowned = *dmmoab->vlocal;
291: #endif
293: /* compute and cache the sizes of local and ghosted entities */
294: dmmoab->nloc = dmmoab->vowned->size();
295: dmmoab->nghost = dmmoab->vghost->size();
296: #ifdef MOAB_HAVE_MPI
297: PetscErrorCode ierr;
298: MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
299: #else
300: dmmoab->n = dmmoab->nloc;
301: #endif
302: return(0);
303: }
306: /*@
307: DMMoabGetAllVertices - Get the entities having DOFs on this DMMoab
309: Collective on MPI_Comm
311: Input Parameter:
312: . dm - The DMMoab object being set
314: Output Parameter:
315: . owned - The local vertex entities in this DMMoab = (owned+ghosted)
317: Level: beginner
319: .keywords: DMMoab, create
320: @*/
321: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
322: {
325: if (local) *local = *((DM_Moab*)dm->data)->vlocal;
326: return(0);
327: }
331: /*@
332: DMMoabGetLocalVertices - Get the entities having DOFs on this DMMoab
334: Collective on MPI_Comm
336: Input Parameter:
337: . dm - The DMMoab object being set
339: Output Parameter:
340: . owned - The owned vertex entities in this DMMoab
341: . ghost - The ghosted entities (non-owned) stored locally in this partition
343: Level: beginner
345: .keywords: DMMoab, create
346: @*/
347: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
348: {
351: if (owned) *owned = ((DM_Moab*)dm->data)->vowned;
352: if (ghost) *ghost = ((DM_Moab*)dm->data)->vghost;
353: return(0);
354: }
356: /*@
357: DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned
359: Collective on MPI_Comm
361: Input Parameter:
362: . dm - The DMMoab object being set
364: Output Parameter:
365: . range - The entities owned locally
367: Level: beginner
369: .keywords: DMMoab, create
370: @*/
371: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
372: {
375: if (range) *range = ((DM_Moab*)dm->data)->elocal;
376: return(0);
377: }
380: /*@
381: DMMoabSetLocalElements - Set the entities having DOFs on this DMMoab
383: Collective on MPI_Comm
385: Input Parameter:
386: . dm - The DMMoab object being set
387: . range - The entities treated by this DMMoab
389: Level: beginner
391: .keywords: DMMoab, create
392: @*/
393: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
394: {
395: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
399: dmmoab->elocal->clear();
400: dmmoab->eghost->clear();
401: dmmoab->elocal->insert(range->begin(), range->end());
402: #ifdef MOAB_HAVE_MPI
403: moab::ErrorCode merr;
404: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
405: *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
406: #endif
407: dmmoab->neleloc = dmmoab->elocal->size();
408: dmmoab->neleghost = dmmoab->eghost->size();
409: #ifdef MOAB_HAVE_MPI
410: PetscErrorCode ierr;
411: MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
412: PetscInfo2(dm, "Created %D local and %D global elements.\n", dmmoab->neleloc, dmmoab->nele);
413: #else
414: dmmoab->nele = dmmoab->neleloc;
415: #endif
416: return(0);
417: }
420: /*@
421: DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
423: Collective on MPI_Comm
425: Input Parameter:
426: . dm - The DMMoab object being set
427: . ltogtag - The MOAB tag used for local to global ids
429: Level: beginner
431: .keywords: DMMoab, create
432: @*/
433: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
434: {
437: ((DM_Moab*)dm->data)->ltog_tag = ltogtag;
438: return(0);
439: }
442: /*@
443: DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
445: Collective on MPI_Comm
447: Input Parameter:
448: . dm - The DMMoab object being set
450: Output Parameter:
451: . ltogtag - The MOAB tag used for local to global ids
453: Level: beginner
455: .keywords: DMMoab, create
456: @*/
457: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
458: {
461: *ltog_tag = ((DM_Moab*)dm->data)->ltog_tag;
462: return(0);
463: }
466: /*@
467: DMMoabSetBlockSize - Set the block size used with this DMMoab
469: Collective on MPI_Comm
471: Input Parameter:
472: . dm - The DMMoab object being set
473: . bs - The block size used with this DMMoab
475: Level: beginner
477: .keywords: DMMoab, create
478: @*/
479: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
480: {
483: ((DM_Moab*)dm->data)->bs = bs;
484: return(0);
485: }
488: /*@
489: DMMoabGetBlockSize - Get the block size used with this DMMoab
491: Collective on MPI_Comm
493: Input Parameter:
494: . dm - The DMMoab object being set
496: Output Parameter:
497: . bs - The block size used with this DMMoab
499: Level: beginner
501: .keywords: DMMoab, create
502: @*/
503: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
504: {
507: *bs = ((DM_Moab*)dm->data)->bs;
508: return(0);
509: }
512: /*@
513: DMMoabGetSize - Get the global vertex size used with this DMMoab
515: Collective on DM
517: Input Parameter:
518: . dm - The DMMoab object being set
520: Output Parameter:
521: . neg - The number of global elements in the DMMoab instance
522: . nvg - The number of global vertices in the DMMoab instance
524: Level: beginner
526: .keywords: DMMoab, create
527: @*/
528: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
529: {
532: if (neg) *neg = ((DM_Moab*)dm->data)->nele;
533: if (nvg) *nvg = ((DM_Moab*)dm->data)->n;
534: return(0);
535: }
538: /*@
539: DMMoabGetLocalSize - Get the local and ghosted vertex size used with this DMMoab
541: Collective on DM
543: Input Parameter:
544: . dm - The DMMoab object being set
546: Output Parameter:
547: + nel - The number of owned elements in this processor
548: . neg - The number of ghosted elements in this processor
549: . nvl - The number of owned vertices in this processor
550: . nvg - The number of ghosted vertices in this processor
552: Level: beginner
554: .keywords: DMMoab, create
555: @*/
556: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
557: {
560: if (nel) *nel = ((DM_Moab*)dm->data)->neleloc;
561: if (neg) *neg = ((DM_Moab*)dm->data)->neleghost;
562: if (nvl) *nvl = ((DM_Moab*)dm->data)->nloc;
563: if (nvg) *nvg = ((DM_Moab*)dm->data)->nghost;
564: return(0);
565: }
568: /*@
569: DMMoabGetOffset - Get the local offset for the global vector
571: Collective on MPI_Comm
573: Input Parameter:
574: . dm - The DMMoab object being set
576: Output Parameter:
577: . offset - The local offset for the global vector
579: Level: beginner
581: .keywords: DMMoab, create
582: @*/
583: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
584: {
587: *offset = ((DM_Moab*)dm->data)->vstart;
588: return(0);
589: }
592: /*@
593: DMMoabGetDimension - Get the dimension of the DM Mesh
595: Collective on MPI_Comm
597: Input Parameter:
598: . dm - The DMMoab object
600: Output Parameter:
601: . dim - The dimension of DM
603: Level: beginner
605: .keywords: DMMoab, create
606: @*/
607: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
608: {
611: *dim = ((DM_Moab*)dm->data)->dim;
612: return(0);
613: }
616: /*@
617: DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
618: generated through uniform refinement.
620: Collective on DM
622: Input Parameter:
623: . dm - The DMMoab object being set
625: Output Parameter:
626: . nvg - The current mesh hierarchy level
628: Level: beginner
630: .keywords: DMMoab, multigrid
631: @*/
632: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
633: {
636: if (nlevel) *nlevel = ((DM_Moab*)dm->data)->hlevel;
637: return(0);
638: }
641: /*@
642: DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh
644: Collective on MPI_Comm
646: Input Parameter:
647: . dm - The DMMoab object
648: . ehandle - The element entity handle
650: Output Parameter:
651: . mat - The material ID for the current entity
653: Level: beginner
655: .keywords: DMMoab, create
656: @*/
657: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
658: {
659: DM_Moab *dmmoab;
663: if (*mat) {
664: dmmoab = (DM_Moab*)(dm)->data;
665: *mat = dmmoab->materials[dmmoab->elocal->index(ehandle)];
666: }
667: return(0);
668: }
671: /*@
672: DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities
674: Collective on MPI_Comm
676: Input Parameter:
677: . dm - The DMMoab object
678: . nconn - Number of entities whose coordinates are needed
679: . conn - The vertex entity handles
681: Output Parameter:
682: . vpos - The coordinates of the requested vertex entities
684: Level: beginner
686: .seealso: DMMoabGetVertexConnectivity()
687: @*/
688: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
689: {
690: DM_Moab *dmmoab;
691: moab::ErrorCode merr;
697: dmmoab = (DM_Moab*)(dm)->data;
699: /* Get connectivity information in MOAB canonical ordering */
700: if (dmmoab->hlevel) {
701: merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle*>(conn), nconn, dmmoab->hlevel, vpos);MBERRNM(merr);
702: }
703: else {
704: merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);MBERRNM(merr);
705: }
706: return(0);
707: }
710: /*@
711: DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity
713: Collective on MPI_Comm
715: Input Parameter:
716: . dm - The DMMoab object
717: . vhandle - Vertex entity handle
719: Output Parameter:
720: . nconn - Number of entities whose coordinates are needed
721: . conn - The vertex entity handles
723: Level: beginner
725: .seealso: DMMoabGetVertexCoordinates(), DMMoabRestoreVertexConnectivity()
726: @*/
727: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt* nconn, moab::EntityHandle **conn)
728: {
729: DM_Moab *dmmoab;
730: std::vector<moab::EntityHandle> adj_entities, connect;
731: PetscErrorCode ierr;
732: moab::ErrorCode merr;
737: dmmoab = (DM_Moab*)(dm)->data;
739: /* Get connectivity information in MOAB canonical ordering */
740: merr = dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION); MBERRNM(merr);
741: merr = dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect); MBERRNM(merr);
743: if (conn) {
744: PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn);
745: PetscMemcpy(*conn, &connect[0], sizeof(moab::EntityHandle) * connect.size());
746: }
747: if (nconn) *nconn = connect.size();
748: return(0);
749: }
752: /*@
753: DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity
755: Collective on MPI_Comm
757: Input Parameter:
758: . dm - The DMMoab object
759: . vhandle - Vertex entity handle
760: . nconn - Number of entities whose coordinates are needed
761: . conn - The vertex entity handles
763: Level: beginner
765: .seealso: DMMoabGetVertexCoordinates(), DMMoabGetVertexConnectivity()
766: @*/
767: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt* nconn, moab::EntityHandle **conn)
768: {
769: PetscErrorCode ierr;
775: if (conn) {
776: PetscFree(*conn);
777: }
778: if (nconn) *nconn = 0;
779: return(0);
780: }
783: /*@
784: DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity
786: Collective on MPI_Comm
788: Input Parameter:
789: . dm - The DMMoab object
790: . ehandle - Vertex entity handle
792: Output Parameter:
793: . nconn - Number of entities whose coordinates are needed
794: . conn - The vertex entity handles
796: Level: beginner
798: .seealso: DMMoabGetVertexCoordinates(), DMMoabGetVertexConnectivity(), DMMoabRestoreVertexConnectivity()
799: @*/
800: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt* nconn, const moab::EntityHandle **conn)
801: {
802: DM_Moab *dmmoab;
803: const moab::EntityHandle *connect;
804: std::vector<moab::EntityHandle> vconn;
805: moab::ErrorCode merr;
806: PetscInt nnodes;
811: dmmoab = (DM_Moab*)(dm)->data;
813: /* Get connectivity information in MOAB canonical ordering */
814: merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes); MBERRNM(merr);
815: if (conn) *conn = connect;
816: if (nconn) *nconn = nnodes;
817: return(0);
818: }
821: /*@
822: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
824: Collective on MPI_Comm
826: Input Parameter:
827: . dm - The DMMoab object
828: . ent - Entity handle
830: Output Parameter:
831: . ent_on_boundary - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
833: Level: beginner
835: .seealso: DMMoabCheckBoundaryVertices()
836: @*/
837: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool* ent_on_boundary)
838: {
839: moab::EntityType etype;
840: DM_Moab *dmmoab;
841: PetscInt edim;
846: dmmoab = (DM_Moab*)(dm)->data;
848: /* get the entity type and handle accordingly */
849: etype = dmmoab->mbiface->type_from_handle(ent);
850: if (etype >= moab::MBPOLYHEDRON) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Entity type on the boundary skin is invalid. EntityType = %D\n", etype);
852: /* get the entity dimension */
853: edim = dmmoab->mbiface->dimension_from_handle(ent);
855: *ent_on_boundary = PETSC_FALSE;
856: if (etype == moab::MBVERTEX && edim == 0) {
857: *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
858: }
859: else {
860: if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
861: if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
862: }
863: else { /* next check the lower-dimensional faces */
864: if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
865: }
866: }
867: return(0);
868: }
871: /*@
872: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
874: Input Parameter:
875: . dm - The DMMoab object
876: . nconn - Number of handles
877: . cnt - Array of entity handles
879: Output Parameter:
880: . isbdvtx - Array of boundary markers - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
882: Level: beginner
884: .seealso: DMMoabIsEntityOnBoundary()
885: @*/
886: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool* isbdvtx)
887: {
888: DM_Moab *dmmoab;
889: PetscInt i;
895: dmmoab = (DM_Moab*)(dm)->data;
897: for (i = 0; i < nconn; ++i) {
898: isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
899: }
900: return(0);
901: }
904: /*@
905: DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary
907: Input Parameter:
908: . dm - The DMMoab object
910: Output Parameter:
911: . bdvtx - Boundary vertices
912: . bdelems - Boundary elements
913: . bdfaces - Boundary faces
915: Level: beginner
917: .seealso: DMMoabCheckBoundaryVertices(), DMMoabIsEntityOnBoundary()
918: @*/
919: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range** bdelems, const moab::Range** bdfaces)
920: {
921: DM_Moab *dmmoab;
925: dmmoab = (DM_Moab*)(dm)->data;
927: if (bdvtx) *bdvtx = dmmoab->bndyvtx;
928: if (bdfaces) *bdfaces = dmmoab->bndyfaces;
929: if (bdelems) *bdfaces = dmmoab->bndyelems;
930: return(0);
931: }
934: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
935: {
936: PetscErrorCode ierr;
937: PetscInt i;
938: moab::ErrorCode merr;
939: DM_Moab *dmmoab = (DM_Moab*)dm->data;
944: dmmoab->refct--;
945: if (!dmmoab->refct) {
946: delete dmmoab->vlocal;
947: delete dmmoab->vowned;
948: delete dmmoab->vghost;
949: delete dmmoab->elocal;
950: delete dmmoab->eghost;
951: delete dmmoab->bndyvtx;
952: delete dmmoab->bndyfaces;
953: delete dmmoab->bndyelems;
955: PetscFree(dmmoab->gsindices);
956: PetscFree2(dmmoab->gidmap, dmmoab->lidmap);
957: PetscFree(dmmoab->dfill);
958: PetscFree(dmmoab->ofill);
959: PetscFree(dmmoab->materials);
960: if (dmmoab->fieldNames) {
961: for (i = 0; i < dmmoab->numFields; i++) {
962: PetscFree(dmmoab->fieldNames[i]);
963: }
964: PetscFree(dmmoab->fieldNames);
965: }
967: if (dmmoab->nhlevels) {
968: PetscFree(dmmoab->hsets);
969: dmmoab->nhlevels = 0;
970: if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
971: dmmoab->hierarchy = NULL;
972: }
974: if (dmmoab->icreatedinstance) {
975: delete dmmoab->pcomm;
976: merr = dmmoab->mbiface->delete_mesh(); MBERRNM(merr);
977: delete dmmoab->mbiface;
978: }
979: dmmoab->mbiface = NULL;
980: #ifdef MOAB_HAVE_MPI
981: dmmoab->pcomm = NULL;
982: #endif
983: VecScatterDestroy(&dmmoab->ltog_sendrecv);
984: ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map);
985: PetscFree(dm->data);
986: }
987: return(0);
988: }
991: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(PetscOptionItems *PetscOptionsObject, DM dm)
992: {
994: DM_Moab *dmmoab = (DM_Moab*)dm->data;
998: PetscOptionsHead(PetscOptionsObject, "DMMoab Options");
999: PetscOptionsInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL);
1000: 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);
1001: /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
1002: 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);
1003: 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);
1004: PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum*)&dmmoab->read_mode, NULL);
1005: PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum*)&dmmoab->write_mode, NULL);
1006: return(0);
1007: }
1010: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
1011: {
1012: PetscErrorCode ierr;
1013: moab::ErrorCode merr;
1014: Vec local, global;
1015: IS from, to;
1016: moab::Range::iterator iter;
1017: PetscInt i, j, f, bs, vent, totsize, *lgmap;
1018: DM_Moab *dmmoab = (DM_Moab*)dm->data;
1019: moab::Range adjs;
1023: /* Get the local and shared vertices and cache it */
1024: if (dmmoab->mbiface == NULL) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB Interface before calling SetUp.");
1025: #ifdef MOAB_HAVE_MPI
1026: if (dmmoab->pcomm == NULL) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB ParallelComm object before calling SetUp.");
1027: #endif
1029: /* Get the entities recursively in the current part of the mesh, if user did not set the local vertices explicitly */
1030: if (dmmoab->vlocal->empty())
1031: {
1032: //merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,*dmmoab->vlocal,true);MBERRNM(merr);
1033: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, 0, *dmmoab->vlocal, false); MBERRNM(merr);
1035: #ifdef MOAB_HAVE_MPI
1036: /* filter based on parallel status */
1037: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned); MBERRNM(merr);
1039: /* filter all the non-owned and shared entities out of the list */
1040: // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
1041: adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
1042: merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost); MBERRNM(merr);
1043: adjs = moab::subtract(adjs, *dmmoab->vghost);
1044: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
1045: #else
1046: *dmmoab->vowned = *dmmoab->vlocal;
1047: #endif
1049: /* compute and cache the sizes of local and ghosted entities */
1050: dmmoab->nloc = dmmoab->vowned->size();
1051: dmmoab->nghost = dmmoab->vghost->size();
1053: #ifdef MOAB_HAVE_MPI
1054: MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
1055: PetscInfo4(NULL, "Filset ID: %u, Vertices: local - %D, owned - %D, ghosted - %D.\n", dmmoab->fileset, dmmoab->vlocal->size(), dmmoab->nloc, dmmoab->nghost);
1056: #else
1057: dmmoab->n = dmmoab->nloc;
1058: #endif
1059: }
1061: {
1062: /* get the information about the local elements in the mesh */
1063: dmmoab->eghost->clear();
1065: /* first decipher the leading dimension */
1066: for (i = 3; i > 0; i--) {
1067: dmmoab->elocal->clear();
1068: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false); MBERRNM(merr);
1070: /* store the current mesh dimension */
1071: if (dmmoab->elocal->size()) {
1072: dmmoab->dim = i;
1073: break;
1074: }
1075: }
1077: DMSetDimension(dm, dmmoab->dim);
1079: #ifdef MOAB_HAVE_MPI
1080: /* filter the ghosted and owned element list */
1081: *dmmoab->eghost = *dmmoab->elocal;
1082: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1083: *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
1084: #endif
1086: dmmoab->neleloc = dmmoab->elocal->size();
1087: dmmoab->neleghost = dmmoab->eghost->size();
1089: #ifdef MOAB_HAVE_MPI
1090: MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
1091: PetscInfo3(NULL, "%d-dim elements: owned - %D, ghosted - %D.\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost);
1092: #else
1093: dmmoab->nele = dmmoab->neleloc;
1094: #endif
1095: }
1097: bs = dmmoab->bs;
1098: if (!dmmoab->ltog_tag) {
1099: /* Get the global ID tag. The global ID tag is applied to each
1100: vertex. It acts as an global identifier which MOAB uses to
1101: assemble the individual pieces of the mesh */
1102: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag); MBERRNM(merr);
1103: }
1105: totsize = dmmoab->vlocal->size();
1106: 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);
1107: PetscCalloc1(totsize, &dmmoab->gsindices);
1108: {
1109: /* first get the local indices */
1110: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]); MBERRNM(merr);
1111: if (dmmoab->nghost) { /* next get the ghosted indices */
1112: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]); MBERRNM(merr);
1113: }
1115: /* find out the local and global minima of GLOBAL_ID */
1116: dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1117: for (i = 0; i < totsize; ++i) {
1118: if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1119: if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1120: }
1122: MPIU_Allreduce(&dmmoab->lminmax[0], &dmmoab->gminmax[0], 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm);
1123: MPIU_Allreduce(&dmmoab->lminmax[1], &dmmoab->gminmax[1], 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm);
1125: /* set the GID map */
1126: for (i = 0; i < totsize; ++i) {
1127: dmmoab->gsindices[i] -= dmmoab->gminmax[0]; /* zero based index needed for IS */
1129: }
1130: dmmoab->lminmax[0] -= dmmoab->gminmax[0];
1131: dmmoab->lminmax[1] -= dmmoab->gminmax[0];
1133: 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]);
1134: }
1135: 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);
1137: {
1138: dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1139: dmmoab->seqend = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1140: PetscInfo2(NULL, "SEQUENCE: Local [min, max] - [%D, %D]\n", dmmoab->seqstart, dmmoab->seqend);
1142: PetscMalloc2(dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->gidmap, dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->lidmap);
1143: PetscMalloc1(totsize * dmmoab->numFields, &lgmap);
1145: i = j = 0;
1146: /* set the owned vertex data first */
1147: for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1148: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1149: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1150: dmmoab->lidmap[vent] = i;
1151: for (f = 0; f < dmmoab->numFields; f++, j++) {
1152: lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1153: }
1154: }
1155: /* next arrange all the ghosted data information */
1156: for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1157: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1158: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1159: dmmoab->lidmap[vent] = i;
1160: for (f = 0; f < dmmoab->numFields; f++, j++) {
1161: lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1162: }
1163: }
1165: /* We need to create the Global to Local Vector Scatter Contexts
1166: 1) First create a local and global vector
1167: 2) Create a local and global IS
1168: 3) Create VecScatter and LtoGMapping objects
1169: 4) Cleanup the IS and Vec objects
1170: */
1171: DMCreateGlobalVector(dm, &global);
1172: DMCreateLocalVector(dm, &local);
1174: VecGetOwnershipRange(global, &dmmoab->vstart, &dmmoab->vend);
1176: /* global to local must retrieve ghost points */
1177: ISCreateStride(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, dmmoab->vstart, 1, &from);
1178: ISSetBlockSize(from, bs);
1180: ISCreateGeneral(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, &lgmap[0], PETSC_COPY_VALUES, &to);
1181: ISSetBlockSize(to, bs);
1183: if (!dmmoab->ltog_map) {
1184: /* create to the local to global mapping for vectors in order to use VecSetValuesLocal */
1185: ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm, dmmoab->bs, totsize * dmmoab->numFields, lgmap,
1186: PETSC_COPY_VALUES, &dmmoab->ltog_map);
1187: }
1189: /* now create the scatter object from local to global vector */
1190: VecScatterCreate(local, from, global, to, &dmmoab->ltog_sendrecv);
1192: /* clean up IS, Vec */
1193: PetscFree(lgmap);
1194: ISDestroy(&from);
1195: ISDestroy(&to);
1196: VecDestroy(&local);
1197: VecDestroy(&global);
1198: }
1200: dmmoab->bndyvtx = new moab::Range();
1201: dmmoab->bndyfaces = new moab::Range();
1202: dmmoab->bndyelems = new moab::Range();
1203: /* skin the boundary and store nodes */
1204: if (!dmmoab->hlevel) {
1205: /* get the skin vertices of boundary faces for the current partition and then filter
1206: the local, boundary faces, vertices and elements alone via PSTATUS flags;
1207: this should not give us any ghosted boundary, but if user needs such a functionality
1208: it would be easy to add it based on the find_skin query below */
1209: moab::Skinner skinner(dmmoab->mbiface);
1211: /* get the entities on the skin - only the faces */
1212: 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
1214: #ifdef MOAB_HAVE_MPI
1215: /* filter all the non-owned and shared entities out of the list */
1216: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1217: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT); MBERRNM(merr);
1218: #endif
1220: /* get all the nodes via connectivity and the parent elements via adjacency information */
1221: merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false); MBERRNM(ierr);
1222: merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION); MBERRNM(ierr);
1223: }
1224: else {
1225: /* Let us query the hierarchy manager and get the results directly for this level */
1226: for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1227: moab::EntityHandle elemHandle = *iter;
1228: if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1229: dmmoab->bndyelems->insert(elemHandle);
1230: /* For this boundary element, query the vertices and add them to the list */
1231: std::vector<moab::EntityHandle> connect;
1232: merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect); MBERRNM(ierr);
1233: for (unsigned iv=0; iv < connect.size(); ++iv)
1234: if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv]))
1235: dmmoab->bndyvtx->insert(connect[iv]);
1236: /* Next, let us query the boundary faces and add them also to the list */
1237: std::vector<moab::EntityHandle> faces;
1238: merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim-1, faces); MBERRNM(ierr);
1239: for (unsigned ifa=0; ifa < faces.size(); ++ifa)
1240: if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa]))
1241: dmmoab->bndyfaces->insert(faces[ifa]);
1242: }
1243: }
1244: #ifdef MOAB_HAVE_MPI
1245: /* filter all the non-owned and shared entities out of the list */
1246: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1247: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1248: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1249: #endif
1251: }
1252: PetscInfo3(NULL, "Found %D boundary vertices, %D boundary faces and %D boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size());
1254: /* Get the material sets and populate the data for all locally owned elements */
1255: {
1256: PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials);
1257: /* Get the count of entities of particular type from dmmoab->elocal
1258: -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1259: moab::Range msets;
1260: 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);
1261: if (msets.size() == 0) {
1262: PetscInfo(NULL, "No material sets found in the fileset.");
1263: }
1265: for (unsigned i=0; i < msets.size(); ++i) {
1266: moab::Range msetelems;
1267: merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);MB_CHK_ERR(merr);
1268: #ifdef MOAB_HAVE_MPI
1269: /* filter all the non-owned and shared entities out of the list */
1270: merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1271: #endif
1273: int partID;
1274: moab::EntityHandle mset=msets[i];
1275: merr = dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID);MB_CHK_ERR(merr);
1277: for (unsigned j=0; j < msetelems.size(); ++j)
1278: dmmoab->materials[dmmoab->elocal->index(msetelems[j])]=partID;
1279: }
1280: }
1282: return(0);
1283: }
1286: /*@
1287: DMMoabCreateVertices - Creates and adds several vertices to the primary set represented by the DM.
1289: Collective on MPI_Comm
1291: Input Parameters:
1292: + dm - The DM object
1293: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1294: . conn - The connectivity of the element
1295: . nverts - The number of vertices that form the element
1297: Output Parameter:
1298: . overts - The list of vertices that were created (can be NULL)
1300: Level: beginner
1302: .keywords: DM, create vertices
1304: .seealso: DMMoabCreateSubmesh(), DMMoabCreateElement()
1305: @*/
1306: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal* coords, PetscInt nverts, moab::Range* overts)
1307: {
1308: moab::ErrorCode merr;
1309: DM_Moab *dmmoab;
1310: moab::Range verts;
1316: dmmoab = (DM_Moab*) dm->data;
1318: /* Insert new points */
1319: merr = dmmoab->mbiface->create_vertices(&coords[0], nverts, verts); MBERRNM(merr);
1320: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, verts); MBERRNM(merr);
1322: if (overts) *overts = verts;
1323: return(0);
1324: }
1327: /*@
1328: DMMoabCreateElement - Adds an element of specified type to the primary set represented by the DM.
1330: Collective on MPI_Comm
1332: Input Parameters:
1333: + dm - The DM object
1334: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1335: . conn - The connectivity of the element
1336: . nverts - The number of vertices that form the element
1338: Output Parameter:
1339: . oelem - The handle to the element created and added to the DM object
1341: Level: beginner
1343: .keywords: DM, create element
1345: .seealso: DMMoabCreateSubmesh(), DMMoabCreateVertices()
1346: @*/
1347: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle* conn, PetscInt nverts, moab::EntityHandle* oelem)
1348: {
1349: moab::ErrorCode merr;
1350: DM_Moab *dmmoab;
1351: moab::EntityHandle elem;
1357: dmmoab = (DM_Moab*) dm->data;
1359: /* Insert new element */
1360: merr = dmmoab->mbiface->create_element(type, conn, nverts, elem); MBERRNM(merr);
1361: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1); MBERRNM(merr);
1363: if (oelem) *oelem = elem;
1364: return(0);
1365: }
1368: /*@
1369: DMMoabCreateSubmesh - Creates a sub-DM object with a set that contains all vertices/elements of the parent
1370: in addition to providing support for dynamic mesh modifications. This is useful for AMR calculations to
1371: create a DM object on a refined level.
1373: Collective on MPI_Comm
1375: Input Parameters:
1376: + dm - The DM object
1378: Output Parameter:
1379: . newdm - The sub DM object with updated set information
1381: Level: advanced
1383: .keywords: DM, sub-DM
1385: .seealso: DMCreate(), DMMoabCreateVertices(), DMMoabCreateElement()
1386: @*/
1387: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1388: {
1389: DM_Moab *dmmoab;
1390: DM_Moab *ndmmoab;
1391: moab::ErrorCode merr;
1392: PetscErrorCode ierr;
1397: dmmoab = (DM_Moab*) dm->data;
1399: /* Create the basic DMMoab object and keep the default parameters created by DM impls */
1400: DMMoabCreateMoab(((PetscObject)dm)->comm, dmmoab->mbiface, &dmmoab->ltog_tag, PETSC_NULL, newdm);
1402: /* get all the necessary handles from the private DM object */
1403: ndmmoab = (DM_Moab*) (*newdm)->data;
1405: /* set the sub-mesh's parent DM reference */
1406: ndmmoab->parent = &dm;
1408: /* create a file set to associate all entities in current mesh */
1409: merr = ndmmoab->mbiface->create_meshset(moab::MESHSET_SET, ndmmoab->fileset); MBERR("Creating file set failed", merr);
1411: /* create a meshset and then add old fileset as child */
1412: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal); MBERR("Adding child vertices to parent failed", merr);
1413: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal); MBERR("Adding child elements to parent failed", merr);
1415: /* preserve the field association between the parent and sub-mesh objects */
1416: DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames);
1417: return(0);
1418: }
1421: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1422: {
1423: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
1424: const char *name;
1425: MPI_Comm comm;
1426: PetscMPIInt size;
1427: PetscErrorCode ierr;
1430: PetscObjectGetComm((PetscObject)dm, &comm);
1431: MPI_Comm_size(comm, &size);
1432: PetscObjectGetName((PetscObject) dm, &name);
1433: PetscViewerASCIIPushTab(viewer);
1434: if (name) {PetscViewerASCIIPrintf(viewer, "%s in %D dimensions:\n", name, dmmoab->dim);}
1435: else {PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dmmoab->dim);}
1436: /* print details about the global mesh */
1437: {
1438: PetscViewerASCIIPushTab(viewer);
1439: PetscViewerASCIIPrintf(viewer, "Sizes: cells=%D, vertices=%D, blocks=%D\n", dmmoab->nele, dmmoab->n, dmmoab->bs);
1440: /* print boundary data */
1441: PetscViewerASCIIPrintf(viewer, "Boundary trace:\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size());
1442: {
1443: PetscViewerASCIIPushTab(viewer);
1444: PetscViewerASCIIPrintf(viewer, "cells=%D, faces=%D, vertices=%D\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size());
1445: PetscViewerASCIIPopTab(viewer);
1446: }
1447: /* print field data */
1448: PetscViewerASCIIPrintf(viewer, "Fields: %D components\n", dmmoab->numFields);
1449: {
1450: PetscViewerASCIIPushTab(viewer);
1451: for (int i = 0; i < dmmoab->numFields; ++i) {
1452: PetscViewerASCIIPrintf(viewer, "[%D] - %s\n", i, dmmoab->fieldNames[i]);
1453: }
1454: PetscViewerASCIIPopTab(viewer);
1455: }
1456: PetscViewerASCIIPopTab(viewer);
1457: }
1458: PetscViewerASCIIPopTab(viewer);
1459: PetscViewerFlush(viewer);
1460: return(0);
1461: }
1463: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1464: {
1465: return(0);
1466: }
1468: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1469: {
1470: return(0);
1471: }
1473: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1474: {
1475: PetscBool iascii, ishdf5, isvtk;
1481: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1482: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
1483: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
1484: if (iascii) {
1485: DMMoabView_Ascii(dm, viewer);
1486: } else if (ishdf5) {
1487: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1488: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
1489: DMMoabView_HDF5(dm, viewer);
1490: PetscViewerPopFormat(viewer);
1491: #else
1492: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1493: #endif
1494: }
1495: else if (isvtk) {
1496: DMMoabView_VTK(dm, viewer);
1497: }
1498: return(0);
1499: }
1502: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
1503: {
1505: dm->ops->view = DMView_Moab;
1506: dm->ops->load = NULL /* DMLoad_Moab */;
1507: dm->ops->setfromoptions = DMSetFromOptions_Moab;
1508: dm->ops->clone = DMClone_Moab;
1509: dm->ops->setup = DMSetUp_Moab;
1510: dm->ops->createdefaultsection = NULL;
1511: dm->ops->createdefaultconstraints = NULL;
1512: dm->ops->createglobalvector = DMCreateGlobalVector_Moab;
1513: dm->ops->createlocalvector = DMCreateLocalVector_Moab;
1514: dm->ops->getlocaltoglobalmapping = NULL;
1515: dm->ops->createfieldis = NULL;
1516: dm->ops->createcoordinatedm = NULL /* DMCreateCoordinateDM_Moab */;
1517: dm->ops->getcoloring = NULL;
1518: dm->ops->creatematrix = DMCreateMatrix_Moab;
1519: dm->ops->createinterpolation = DMCreateInterpolation_Moab;
1520: dm->ops->getaggregates = NULL;
1521: dm->ops->getinjection = NULL /* DMCreateInjection_Moab */;
1522: dm->ops->refine = DMRefine_Moab;
1523: dm->ops->coarsen = DMCoarsen_Moab;
1524: dm->ops->refinehierarchy = DMRefineHierarchy_Moab;
1525: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Moab;
1526: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Moab;
1527: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Moab;
1528: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Moab;
1529: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Moab;
1530: dm->ops->destroy = DMDestroy_Moab;
1531: dm->ops->createsubdm = NULL /* DMCreateSubDM_Moab */;
1532: dm->ops->getdimpoints = NULL /* DMGetDimPoints_Moab */;
1533: dm->ops->locatepoints = NULL /* DMLocatePoints_Moab */;
1534: return(0);
1535: }
1538: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1539: {
1540: PetscErrorCode ierr;
1543: PetscObjectChangeTypeName((PetscObject) * newdm, DMMOAB);
1545: /* get all the necessary handles from the private DM object */
1546: (*newdm)->data = (DM_Moab*) dm->data;
1547: ((DM_Moab*)dm->data)->refct++;
1549: DMInitialize_Moab(*newdm);
1550: return(0);
1551: }
1554: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1555: {
1560: PetscNewLog(dm, (DM_Moab**)&dm->data);
1562: ((DM_Moab*)dm->data)->bs = 1;
1563: ((DM_Moab*)dm->data)->numFields = 1;
1564: ((DM_Moab*)dm->data)->n = 0;
1565: ((DM_Moab*)dm->data)->nloc = 0;
1566: ((DM_Moab*)dm->data)->nghost = 0;
1567: ((DM_Moab*)dm->data)->nele = 0;
1568: ((DM_Moab*)dm->data)->neleloc = 0;
1569: ((DM_Moab*)dm->data)->neleghost = 0;
1570: ((DM_Moab*)dm->data)->ltog_map = NULL;
1571: ((DM_Moab*)dm->data)->ltog_sendrecv = NULL;
1573: ((DM_Moab*)dm->data)->refct = 1;
1574: ((DM_Moab*)dm->data)->parent = NULL;
1575: ((DM_Moab*)dm->data)->vlocal = new moab::Range();
1576: ((DM_Moab*)dm->data)->vowned = new moab::Range();
1577: ((DM_Moab*)dm->data)->vghost = new moab::Range();
1578: ((DM_Moab*)dm->data)->elocal = new moab::Range();
1579: ((DM_Moab*)dm->data)->eghost = new moab::Range();
1581: DMInitialize_Moab(dm);
1582: return(0);
1583: }