Actual source code: dmmoab.cxx
1: #include <petsc/private/dmmbimpl.h>
3: #include <petscdmmoab.h>
4: #include <MBTagConventions.hpp>
5: #include <moab/NestedRefine.hpp>
6: #include <moab/Skinner.hpp>
8: /*MC
9: DMMOAB = "moab" - A DM object that encapsulates an unstructured mesh described by the MOAB mesh database.
10: Direct access to the MOAB Interface and other mesh manipulation related objects are available
11: through public API. Ability to create global and local representation of Vecs containing all
12: unknowns in the interior and shared boundary via a transparent tag-data wrapper is provided
13: along with utility functions to traverse the mesh and assemble a discrete system via
14: field-based/blocked Vec(Get/Set) methods. Input from and output to different formats are
15: available.
17: Reference: https://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html
19: Level: intermediate
21: .seealso: DMType, DMMoabCreate(), DMCreate(), DMSetType(), DMMoabCreateMoab()
22: M*/
24: /* External function declarations here */
25: PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
26: PETSC_EXTERN PetscErrorCode DMCreateDefaultConstraints_Moab(DM dm);
27: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
28: PETSC_EXTERN PetscErrorCode DMCreateCoordinateDM_Moab(DM dm, DM *cdm);
29: PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM *dmRefined);
30: PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM *dmCoarsened);
31: PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmRefined[]);
32: PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmCoarsened[]);
33: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm);
34: PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM, Vec *);
35: PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM, Vec *);
36: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
37: PETSC_EXTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM, Vec, InsertMode, Vec);
38: PETSC_EXTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM, Vec, InsertMode, Vec);
39: PETSC_EXTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM, Vec, InsertMode, Vec);
40: PETSC_EXTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM, Vec, InsertMode, Vec);
42: /* Un-implemented routines */
43: /*
44: PETSC_EXTERN PetscErrorCode DMCreatelocalsection_Moab(DM dm);
45: PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dmCoarse, DM dmFine, Mat *mat);
46: PETSC_EXTERN PetscErrorCode DMLoad_Moab(DM dm, PetscViewer viewer);
47: PETSC_EXTERN PetscErrorCode DMGetDimPoints_Moab(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd);
48: PETSC_EXTERN PetscErrorCode DMCreateSubDM_Moab(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
49: PETSC_EXTERN PetscErrorCode DMLocatePoints_Moab(DM dm, Vec v, IS *cellIS);
50: */
52: /*@C
53: DMMoabCreate - Creates a DMMoab object, which encapsulates a moab instance
55: Collective
57: Input Parameter:
58: . comm - The communicator for the DMMoab object
60: Output Parameter:
61: . dmb - The DMMoab object
63: Level: beginner
65: @*/
66: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
67: {
72: DMCreate(comm, dmb);
73: DMSetType(*dmb, DMMOAB);
74: return(0);
75: }
77: /*@C
78: DMMoabCreateMoab - Creates a DMMoab object, optionally from an instance and other data
80: Collective
82: Input Parameters:
83: + comm - The communicator for the DMMoab object
84: . mbiface - (ptr to) the MOAB Instance; if passed in NULL, MOAB instance is created inside PETSc, and destroyed
85: along with the DMMoab
86: . ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
87: - range - If non-NULL, contains range of entities to which DOFs will be assigned
89: Output Parameter:
90: . dmb - The DMMoab object
92: Level: intermediate
94: @*/
95: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
96: {
98: moab::ErrorCode merr;
99: DM dmmb;
100: DM_Moab *dmmoab;
105: DMMoabCreate(comm, &dmmb);
106: dmmoab = (DM_Moab*)(dmmb)->data;
108: if (!mbiface) {
109: dmmoab->mbiface = new moab::Core();
110: dmmoab->icreatedinstance = PETSC_TRUE;
111: }
112: else {
113: dmmoab->mbiface = mbiface;
114: dmmoab->icreatedinstance = PETSC_FALSE;
115: }
117: /* by default the fileset = root set. This set stores the hierarchy of entities belonging to current DM */
118: dmmoab->fileset = 0;
119: dmmoab->hlevel = 0;
120: dmmoab->nghostrings = 0;
122: #ifdef MOAB_HAVE_MPI
123: moab::EntityHandle partnset;
125: /* Create root sets for each mesh. Then pass these
126: to the load_file functions to be populated. */
127: merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, partnset); MBERR("Creating partition set failed", merr);
129: /* Create the parallel communicator object with the partition handle associated with MOAB */
130: dmmoab->pcomm = moab::ParallelComm::get_pcomm(dmmoab->mbiface, partnset, &comm);
131: #endif
133: /* do the remaining initializations for DMMoab */
134: dmmoab->bs = 1;
135: dmmoab->numFields = 1;
136: PetscMalloc(dmmoab->numFields * sizeof(char*), &dmmoab->fieldNames);
137: PetscStrallocpy("DEFAULT", (char**) &dmmoab->fieldNames[0]);
138: dmmoab->rw_dbglevel = 0;
139: dmmoab->partition_by_rank = PETSC_FALSE;
140: dmmoab->extra_read_options[0] = '\0';
141: dmmoab->extra_write_options[0] = '\0';
142: dmmoab->read_mode = READ_PART;
143: dmmoab->write_mode = WRITE_PART;
145: /* set global ID tag handle */
146: if (ltog_tag && *ltog_tag) {
147: DMMoabSetLocalToGlobalTag(dmmb, *ltog_tag);
148: }
149: else {
150: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag); MBERRNM(merr);
151: if (ltog_tag) *ltog_tag = dmmoab->ltog_tag;
152: }
154: merr = dmmoab->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dmmoab->material_tag); MBERRNM(merr);
156: /* set the local range of entities (vertices) of interest */
157: if (range) {
158: DMMoabSetLocalVertices(dmmb, range);
159: }
160: *dmb = dmmb;
161: return(0);
162: }
164: #ifdef MOAB_HAVE_MPI
166: /*@C
167: DMMoabGetParallelComm - Get the ParallelComm used with this DMMoab
169: Collective
171: Input Parameter:
172: . dm - The DMMoab object being set
174: Output Parameter:
175: . pcomm - The ParallelComm for the DMMoab
177: Level: beginner
179: @*/
180: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
181: {
184: *pcomm = ((DM_Moab*)(dm)->data)->pcomm;
185: return(0);
186: }
188: #endif /* MOAB_HAVE_MPI */
190: /*@C
191: DMMoabSetInterface - Set the MOAB instance used with this DMMoab
193: Collective
195: Input Parameters:
196: + dm - The DMMoab object being set
197: - mbiface - The MOAB instance being set on this DMMoab
199: Level: beginner
201: @*/
202: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
203: {
204: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
209: #ifdef MOAB_HAVE_MPI
210: dmmoab->pcomm = NULL;
211: #endif
212: dmmoab->mbiface = mbiface;
213: dmmoab->icreatedinstance = PETSC_FALSE;
214: return(0);
215: }
217: /*@C
218: DMMoabGetInterface - Get the MOAB instance used with this DMMoab
220: Collective
222: Input Parameter:
223: . dm - The DMMoab object being set
225: Output Parameter:
226: . mbiface - The MOAB instance set on this DMMoab
228: Level: beginner
230: @*/
231: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
232: {
233: PetscErrorCode ierr;
234: static PetscBool cite = PETSC_FALSE;
238: 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);
239: *mbiface = ((DM_Moab*)dm->data)->mbiface;
240: return(0);
241: }
243: /*@C
244: DMMoabSetLocalVertices - Set the entities having DOFs on this DMMoab
246: Collective
248: Input Parameters:
249: + dm - The DMMoab object being set
250: - range - The entities treated by this DMMoab
252: Level: beginner
254: @*/
255: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
256: {
257: moab::Range tmpvtxs;
258: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
262: dmmoab->vlocal->clear();
263: dmmoab->vowned->clear();
265: dmmoab->vlocal->insert(range->begin(), range->end());
267: #ifdef MOAB_HAVE_MPI
268: moab::ErrorCode merr;
269: /* filter based on parallel status */
270: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned); MBERRNM(merr);
272: /* filter all the non-owned and shared entities out of the list */
273: tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
274: merr = dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost); MBERRNM(merr);
275: tmpvtxs = moab::subtract(tmpvtxs, *dmmoab->vghost);
276: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
277: #else
278: *dmmoab->vowned = *dmmoab->vlocal;
279: #endif
281: /* compute and cache the sizes of local and ghosted entities */
282: dmmoab->nloc = dmmoab->vowned->size();
283: dmmoab->nghost = dmmoab->vghost->size();
284: #ifdef MOAB_HAVE_MPI
285: PetscErrorCode ierr;
286: MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
287: #else
288: dmmoab->n = dmmoab->nloc;
289: #endif
290: return(0);
291: }
293: /*@C
294: DMMoabGetAllVertices - Get the entities having DOFs on this DMMoab
296: Collective
298: Input Parameter:
299: . dm - The DMMoab object being set
301: Output Parameter:
302: . owned - The local vertex entities in this DMMoab = (owned+ghosted)
304: Level: beginner
306: @*/
307: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
308: {
311: if (local) *local = *((DM_Moab*)dm->data)->vlocal;
312: return(0);
313: }
315: /*@C
316: DMMoabGetLocalVertices - Get the entities having DOFs on this DMMoab
318: Collective
320: Input Parameter:
321: . dm - The DMMoab object being set
323: Output Parameters:
324: + owned - The owned vertex entities in this DMMoab
325: - ghost - The ghosted entities (non-owned) stored locally in this partition
327: Level: beginner
329: @*/
330: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
331: {
334: if (owned) *owned = ((DM_Moab*)dm->data)->vowned;
335: if (ghost) *ghost = ((DM_Moab*)dm->data)->vghost;
336: return(0);
337: }
339: /*@C
340: DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned
342: Collective
344: Input Parameter:
345: . dm - The DMMoab object being set
347: Output Parameter:
348: . range - The entities owned locally
350: Level: beginner
352: @*/
353: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
354: {
357: if (range) *range = ((DM_Moab*)dm->data)->elocal;
358: return(0);
359: }
361: /*@C
362: DMMoabSetLocalElements - Set the entities having DOFs on this DMMoab
364: Collective
366: Input Parameters:
367: + dm - The DMMoab object being set
368: - range - The entities treated by this DMMoab
370: Level: beginner
372: @*/
373: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
374: {
375: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
379: dmmoab->elocal->clear();
380: dmmoab->eghost->clear();
381: dmmoab->elocal->insert(range->begin(), range->end());
382: #ifdef MOAB_HAVE_MPI
383: moab::ErrorCode merr;
384: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
385: *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
386: #endif
387: dmmoab->neleloc = dmmoab->elocal->size();
388: dmmoab->neleghost = dmmoab->eghost->size();
389: #ifdef MOAB_HAVE_MPI
390: PetscErrorCode ierr;
391: MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
392: PetscInfo2(dm, "Created %D local and %D global elements.\n", dmmoab->neleloc, dmmoab->nele);
393: #else
394: dmmoab->nele = dmmoab->neleloc;
395: #endif
396: return(0);
397: }
399: /*@C
400: DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
402: Collective
404: Input Parameters:
405: + dm - The DMMoab object being set
406: - ltogtag - The MOAB tag used for local to global ids
408: Level: beginner
410: @*/
411: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
412: {
415: ((DM_Moab*)dm->data)->ltog_tag = ltogtag;
416: return(0);
417: }
419: /*@C
420: DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
422: Collective
424: Input Parameter:
425: . dm - The DMMoab object being set
427: Output Parameter:
428: . ltogtag - The MOAB tag used for local to global ids
430: Level: beginner
432: @*/
433: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
434: {
437: *ltog_tag = ((DM_Moab*)dm->data)->ltog_tag;
438: return(0);
439: }
441: /*@C
442: DMMoabSetBlockSize - Set the block size used with this DMMoab
444: Collective
446: Input Parameters:
447: + dm - The DMMoab object being set
448: - bs - The block size used with this DMMoab
450: Level: beginner
452: @*/
453: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
454: {
457: ((DM_Moab*)dm->data)->bs = bs;
458: return(0);
459: }
461: /*@C
462: DMMoabGetBlockSize - Get the block size used with this DMMoab
464: Collective
466: Input Parameter:
467: . dm - The DMMoab object being set
469: Output Parameter:
470: . bs - The block size used with this DMMoab
472: Level: beginner
474: @*/
475: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
476: {
479: *bs = ((DM_Moab*)dm->data)->bs;
480: return(0);
481: }
483: /*@C
484: DMMoabGetSize - Get the global vertex size used with this DMMoab
486: Collective on dm
488: Input Parameter:
489: . dm - The DMMoab object being set
491: Output Parameters:
492: + neg - The number of global elements in the DMMoab instance
493: - nvg - The number of global vertices in the DMMoab instance
495: Level: beginner
497: @*/
498: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
499: {
502: if (neg) *neg = ((DM_Moab*)dm->data)->nele;
503: if (nvg) *nvg = ((DM_Moab*)dm->data)->n;
504: return(0);
505: }
507: /*@C
508: DMMoabGetLocalSize - Get the local and ghosted vertex size used with this DMMoab
510: Collective on dm
512: Input Parameter:
513: . dm - The DMMoab object being set
515: Output Parameters:
516: + nel - The number of owned elements in this processor
517: . neg - The number of ghosted elements in this processor
518: . nvl - The number of owned vertices in this processor
519: - nvg - The number of ghosted vertices in this processor
521: Level: beginner
523: @*/
524: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
525: {
528: if (nel) *nel = ((DM_Moab*)dm->data)->neleloc;
529: if (neg) *neg = ((DM_Moab*)dm->data)->neleghost;
530: if (nvl) *nvl = ((DM_Moab*)dm->data)->nloc;
531: if (nvg) *nvg = ((DM_Moab*)dm->data)->nghost;
532: return(0);
533: }
535: /*@C
536: DMMoabGetOffset - Get the local offset for the global vector
538: Collective
540: Input Parameter:
541: . dm - The DMMoab object being set
543: Output Parameter:
544: . offset - The local offset for the global vector
546: Level: beginner
548: @*/
549: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
550: {
553: *offset = ((DM_Moab*)dm->data)->vstart;
554: return(0);
555: }
557: /*@C
558: DMMoabGetDimension - Get the dimension of the DM Mesh
560: Collective
562: Input Parameter:
563: . dm - The DMMoab object
565: Output Parameter:
566: . dim - The dimension of DM
568: Level: beginner
570: @*/
571: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
572: {
575: *dim = ((DM_Moab*)dm->data)->dim;
576: return(0);
577: }
579: /*@C
580: DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
581: generated through uniform refinement.
583: Collective on dm
585: Input Parameter:
586: . dm - The DMMoab object being set
588: Output Parameter:
589: . nvg - The current mesh hierarchy level
591: Level: beginner
593: @*/
594: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
595: {
598: if (nlevel) *nlevel = ((DM_Moab*)dm->data)->hlevel;
599: return(0);
600: }
602: /*@C
603: DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh
605: Collective
607: Input Parameters:
608: + dm - The DMMoab object
609: - ehandle - The element entity handle
611: Output Parameter:
612: . mat - The material ID for the current entity
614: Level: beginner
616: @*/
617: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
618: {
619: DM_Moab *dmmoab;
623: if (*mat) {
624: dmmoab = (DM_Moab*)(dm)->data;
625: *mat = dmmoab->materials[dmmoab->elocal->index(ehandle)];
626: }
627: return(0);
628: }
630: /*@C
631: DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities
633: Collective
635: Input Parameters:
636: + dm - The DMMoab object
637: . nconn - Number of entities whose coordinates are needed
638: - conn - The vertex entity handles
640: Output Parameter:
641: . vpos - The coordinates of the requested vertex entities
643: Level: beginner
645: .seealso: DMMoabGetVertexConnectivity()
646: @*/
647: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
648: {
649: DM_Moab *dmmoab;
650: moab::ErrorCode merr;
656: dmmoab = (DM_Moab*)(dm)->data;
658: /* Get connectivity information in MOAB canonical ordering */
659: if (dmmoab->hlevel) {
660: merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle*>(conn), nconn, dmmoab->hlevel, vpos);MBERRNM(merr);
661: }
662: else {
663: merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);MBERRNM(merr);
664: }
665: return(0);
666: }
668: /*@C
669: DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity
671: Collective
673: Input Parameters:
674: + dm - The DMMoab object
675: - vhandle - Vertex entity handle
677: Output Parameters:
678: + nconn - Number of entities whose coordinates are needed
679: - conn - The vertex entity handles
681: Level: beginner
683: .seealso: DMMoabGetVertexCoordinates(), DMMoabRestoreVertexConnectivity()
684: @*/
685: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt* nconn, moab::EntityHandle **conn)
686: {
687: DM_Moab *dmmoab;
688: std::vector<moab::EntityHandle> adj_entities, connect;
689: PetscErrorCode ierr;
690: moab::ErrorCode merr;
695: dmmoab = (DM_Moab*)(dm)->data;
697: /* Get connectivity information in MOAB canonical ordering */
698: merr = dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION); MBERRNM(merr);
699: merr = dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect); MBERRNM(merr);
701: if (conn) {
702: PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn);
703: PetscArraycpy(*conn, &connect[0], connect.size());
704: }
705: if (nconn) *nconn = connect.size();
706: return(0);
707: }
709: /*@C
710: DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity
712: Collective
714: Input Parameters:
715: + dm - The DMMoab object
716: . vhandle - Vertex entity handle
717: . nconn - Number of entities whose coordinates are needed
718: - conn - The vertex entity handles
720: Level: beginner
722: .seealso: DMMoabGetVertexCoordinates(), DMMoabGetVertexConnectivity()
723: @*/
724: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt* nconn, moab::EntityHandle **conn)
725: {
726: PetscErrorCode ierr;
732: if (conn) {
733: PetscFree(*conn);
734: }
735: if (nconn) *nconn = 0;
736: return(0);
737: }
739: /*@C
740: DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity
742: Collective
744: Input Parameters:
745: + dm - The DMMoab object
746: - ehandle - Vertex entity handle
748: Output Parameters:
749: + nconn - Number of entities whose coordinates are needed
750: - conn - The vertex entity handles
752: Level: beginner
754: .seealso: DMMoabGetVertexCoordinates(), DMMoabGetVertexConnectivity(), DMMoabRestoreVertexConnectivity()
755: @*/
756: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt* nconn, const moab::EntityHandle **conn)
757: {
758: DM_Moab *dmmoab;
759: const moab::EntityHandle *connect;
760: std::vector<moab::EntityHandle> vconn;
761: moab::ErrorCode merr;
762: PetscInt nnodes;
767: dmmoab = (DM_Moab*)(dm)->data;
769: /* Get connectivity information in MOAB canonical ordering */
770: merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes); MBERRNM(merr);
771: if (conn) *conn = connect;
772: if (nconn) *nconn = nnodes;
773: return(0);
774: }
776: /*@C
777: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
779: Collective
781: Input Parameters:
782: + dm - The DMMoab object
783: - ent - Entity handle
785: Output Parameter:
786: . ent_on_boundary - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
788: Level: beginner
790: .seealso: DMMoabCheckBoundaryVertices()
791: @*/
792: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool* ent_on_boundary)
793: {
794: moab::EntityType etype;
795: DM_Moab *dmmoab;
796: PetscInt edim;
801: dmmoab = (DM_Moab*)(dm)->data;
803: /* get the entity type and handle accordingly */
804: etype = dmmoab->mbiface->type_from_handle(ent);
805: if (etype >= moab::MBPOLYHEDRON) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Entity type on the boundary skin is invalid. EntityType = %D\n", etype);
807: /* get the entity dimension */
808: edim = dmmoab->mbiface->dimension_from_handle(ent);
810: *ent_on_boundary = PETSC_FALSE;
811: if (etype == moab::MBVERTEX && edim == 0) {
812: *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
813: }
814: else {
815: if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
816: if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
817: }
818: else { /* next check the lower-dimensional faces */
819: if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
820: }
821: }
822: return(0);
823: }
825: /*@C
826: DMMoabCheckBoundaryVertices - Check whether a given entity is on the boundary (vertex, edge, face, element)
828: Input Parameters:
829: + dm - The DMMoab object
830: . nconn - Number of handles
831: - cnt - Array of entity handles
833: Output Parameter:
834: . isbdvtx - Array of boundary markers - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
836: Level: beginner
838: .seealso: DMMoabIsEntityOnBoundary()
839: @*/
840: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool* isbdvtx)
841: {
842: DM_Moab *dmmoab;
843: PetscInt i;
849: dmmoab = (DM_Moab*)(dm)->data;
851: for (i = 0; i < nconn; ++i) {
852: isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
853: }
854: return(0);
855: }
857: /*@C
858: DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary
860: Input Parameter:
861: . dm - The DMMoab object
863: Output Parameters:
864: + bdvtx - Boundary vertices
865: . bdelems - Boundary elements
866: - bdfaces - Boundary faces
868: Level: beginner
870: .seealso: DMMoabCheckBoundaryVertices(), DMMoabIsEntityOnBoundary()
871: @*/
872: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range** bdelems, const moab::Range** bdfaces)
873: {
874: DM_Moab *dmmoab;
878: dmmoab = (DM_Moab*)(dm)->data;
880: if (bdvtx) *bdvtx = dmmoab->bndyvtx;
881: if (bdfaces) *bdfaces = dmmoab->bndyfaces;
882: if (bdelems) *bdfaces = dmmoab->bndyelems;
883: return(0);
884: }
886: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
887: {
888: PetscErrorCode ierr;
889: PetscInt i;
890: moab::ErrorCode merr;
891: DM_Moab *dmmoab = (DM_Moab*)dm->data;
896: dmmoab->refct--;
897: if (!dmmoab->refct) {
898: delete dmmoab->vlocal;
899: delete dmmoab->vowned;
900: delete dmmoab->vghost;
901: delete dmmoab->elocal;
902: delete dmmoab->eghost;
903: delete dmmoab->bndyvtx;
904: delete dmmoab->bndyfaces;
905: delete dmmoab->bndyelems;
907: PetscFree(dmmoab->gsindices);
908: PetscFree2(dmmoab->gidmap, dmmoab->lidmap);
909: PetscFree(dmmoab->dfill);
910: PetscFree(dmmoab->ofill);
911: PetscFree(dmmoab->materials);
912: if (dmmoab->fieldNames) {
913: for (i = 0; i < dmmoab->numFields; i++) {
914: PetscFree(dmmoab->fieldNames[i]);
915: }
916: PetscFree(dmmoab->fieldNames);
917: }
919: if (dmmoab->nhlevels) {
920: PetscFree(dmmoab->hsets);
921: dmmoab->nhlevels = 0;
922: if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
923: dmmoab->hierarchy = NULL;
924: }
926: if (dmmoab->icreatedinstance) {
927: delete dmmoab->pcomm;
928: merr = dmmoab->mbiface->delete_mesh(); MBERRNM(merr);
929: delete dmmoab->mbiface;
930: }
931: dmmoab->mbiface = NULL;
932: #ifdef MOAB_HAVE_MPI
933: dmmoab->pcomm = NULL;
934: #endif
935: VecScatterDestroy(&dmmoab->ltog_sendrecv);
936: ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map);
937: PetscFree(dm->data);
938: }
939: return(0);
940: }
942: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(PetscOptionItems *PetscOptionsObject, DM dm)
943: {
945: DM_Moab *dmmoab = (DM_Moab*)dm->data;
949: PetscOptionsHead(PetscOptionsObject, "DMMoab Options");
950: PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL,0);
951: 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);
952: /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
953: PetscOptionsString("-dm_moab_read_opts", "Extra options to enable MOAB reader to load DM from file", "DMView", dmmoab->extra_read_options, dmmoab->extra_read_options, sizeof(dmmoab->extra_read_options), NULL);
954: PetscOptionsString("-dm_moab_write_opts", "Extra options to enable MOAB writer to serialize DM to file", "DMView", dmmoab->extra_write_options, dmmoab->extra_write_options, sizeof(dmmoab->extra_write_options), NULL);
955: PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum*)&dmmoab->read_mode, NULL);
956: PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum*)&dmmoab->write_mode, NULL);
957: return(0);
958: }
960: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
961: {
962: PetscErrorCode ierr;
963: moab::ErrorCode merr;
964: Vec local, global;
965: IS from, to;
966: moab::Range::iterator iter;
967: PetscInt i, j, f, bs, vent, totsize, *lgmap;
968: DM_Moab *dmmoab = (DM_Moab*)dm->data;
969: moab::Range adjs;
973: /* Get the local and shared vertices and cache it */
974: if (dmmoab->mbiface == NULL) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB Interface before calling SetUp.");
975: #ifdef MOAB_HAVE_MPI
976: if (dmmoab->pcomm == NULL) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB ParallelComm object before calling SetUp.");
977: #endif
979: /* Get the entities recursively in the current part of the mesh, if user did not set the local vertices explicitly */
980: if (dmmoab->vlocal->empty())
981: {
982: //merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,*dmmoab->vlocal,true);MBERRNM(merr);
983: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, 0, *dmmoab->vlocal, false); MBERRNM(merr);
985: #ifdef MOAB_HAVE_MPI
986: /* filter based on parallel status */
987: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned); MBERRNM(merr);
989: /* filter all the non-owned and shared entities out of the list */
990: // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
991: adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
992: merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost); MBERRNM(merr);
993: adjs = moab::subtract(adjs, *dmmoab->vghost);
994: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
995: #else
996: *dmmoab->vowned = *dmmoab->vlocal;
997: #endif
999: /* compute and cache the sizes of local and ghosted entities */
1000: dmmoab->nloc = dmmoab->vowned->size();
1001: dmmoab->nghost = dmmoab->vghost->size();
1003: #ifdef MOAB_HAVE_MPI
1004: MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
1005: PetscInfo4(NULL, "Filset ID: %u, Vertices: local - %D, owned - %D, ghosted - %D.\n", dmmoab->fileset, dmmoab->vlocal->size(), dmmoab->nloc, dmmoab->nghost);
1006: #else
1007: dmmoab->n = dmmoab->nloc;
1008: #endif
1009: }
1011: {
1012: /* get the information about the local elements in the mesh */
1013: dmmoab->eghost->clear();
1015: /* first decipher the leading dimension */
1016: for (i = 3; i > 0; i--) {
1017: dmmoab->elocal->clear();
1018: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false); MBERRNM(merr);
1020: /* store the current mesh dimension */
1021: if (dmmoab->elocal->size()) {
1022: dmmoab->dim = i;
1023: break;
1024: }
1025: }
1027: DMSetDimension(dm, dmmoab->dim);
1029: #ifdef MOAB_HAVE_MPI
1030: /* filter the ghosted and owned element list */
1031: *dmmoab->eghost = *dmmoab->elocal;
1032: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1033: *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
1034: #endif
1036: dmmoab->neleloc = dmmoab->elocal->size();
1037: dmmoab->neleghost = dmmoab->eghost->size();
1039: #ifdef MOAB_HAVE_MPI
1040: MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
1041: PetscInfo3(NULL, "%d-dim elements: owned - %D, ghosted - %D.\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost);
1042: #else
1043: dmmoab->nele = dmmoab->neleloc;
1044: #endif
1045: }
1047: bs = dmmoab->bs;
1048: if (!dmmoab->ltog_tag) {
1049: /* Get the global ID tag. The global ID tag is applied to each
1050: vertex. It acts as an global identifier which MOAB uses to
1051: assemble the individual pieces of the mesh */
1052: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag); MBERRNM(merr);
1053: }
1055: totsize = dmmoab->vlocal->size();
1056: 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);
1057: PetscCalloc1(totsize, &dmmoab->gsindices);
1058: {
1059: /* first get the local indices */
1060: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]); MBERRNM(merr);
1061: if (dmmoab->nghost) { /* next get the ghosted indices */
1062: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]); MBERRNM(merr);
1063: }
1065: /* find out the local and global minima of GLOBAL_ID */
1066: dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1067: for (i = 0; i < totsize; ++i) {
1068: if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1069: if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1070: }
1072: MPIU_Allreduce(&dmmoab->lminmax[0], &dmmoab->gminmax[0], 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm);
1073: MPIU_Allreduce(&dmmoab->lminmax[1], &dmmoab->gminmax[1], 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm);
1075: /* set the GID map */
1076: for (i = 0; i < totsize; ++i) {
1077: dmmoab->gsindices[i] -= dmmoab->gminmax[0]; /* zero based index needed for IS */
1079: }
1080: dmmoab->lminmax[0] -= dmmoab->gminmax[0];
1081: dmmoab->lminmax[1] -= dmmoab->gminmax[0];
1083: 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]);
1084: }
1085: 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);
1087: {
1088: dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1089: dmmoab->seqend = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1090: PetscInfo2(NULL, "SEQUENCE: Local [min, max] - [%D, %D]\n", dmmoab->seqstart, dmmoab->seqend);
1092: PetscMalloc2(dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->gidmap, dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->lidmap);
1093: PetscMalloc1(totsize * dmmoab->numFields, &lgmap);
1095: i = j = 0;
1096: /* set the owned vertex data first */
1097: for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1098: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1099: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1100: dmmoab->lidmap[vent] = i;
1101: for (f = 0; f < dmmoab->numFields; f++, j++) {
1102: lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1103: }
1104: }
1105: /* next arrange all the ghosted data information */
1106: for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1107: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1108: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1109: dmmoab->lidmap[vent] = i;
1110: for (f = 0; f < dmmoab->numFields; f++, j++) {
1111: lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1112: }
1113: }
1115: /* We need to create the Global to Local Vector Scatter Contexts
1116: 1) First create a local and global vector
1117: 2) Create a local and global IS
1118: 3) Create VecScatter and LtoGMapping objects
1119: 4) Cleanup the IS and Vec objects
1120: */
1121: DMCreateGlobalVector(dm, &global);
1122: DMCreateLocalVector(dm, &local);
1124: VecGetOwnershipRange(global, &dmmoab->vstart, &dmmoab->vend);
1126: /* global to local must retrieve ghost points */
1127: ISCreateStride(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, dmmoab->vstart, 1, &from);
1128: ISSetBlockSize(from, bs);
1130: ISCreateGeneral(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, &lgmap[0], PETSC_COPY_VALUES, &to);
1131: ISSetBlockSize(to, bs);
1133: if (!dmmoab->ltog_map) {
1134: /* create to the local to global mapping for vectors in order to use VecSetValuesLocal */
1135: ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm, dmmoab->bs, totsize * dmmoab->numFields, lgmap,
1136: PETSC_COPY_VALUES, &dmmoab->ltog_map);
1137: }
1139: /* now create the scatter object from local to global vector */
1140: VecScatterCreate(local, from, global, to, &dmmoab->ltog_sendrecv);
1142: /* clean up IS, Vec */
1143: PetscFree(lgmap);
1144: ISDestroy(&from);
1145: ISDestroy(&to);
1146: VecDestroy(&local);
1147: VecDestroy(&global);
1148: }
1150: dmmoab->bndyvtx = new moab::Range();
1151: dmmoab->bndyfaces = new moab::Range();
1152: dmmoab->bndyelems = new moab::Range();
1153: /* skin the boundary and store nodes */
1154: if (!dmmoab->hlevel) {
1155: /* get the skin vertices of boundary faces for the current partition and then filter
1156: the local, boundary faces, vertices and elements alone via PSTATUS flags;
1157: this should not give us any ghosted boundary, but if user needs such a functionality
1158: it would be easy to add it based on the find_skin query below */
1159: moab::Skinner skinner(dmmoab->mbiface);
1161: /* get the entities on the skin - only the faces */
1162: 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
1164: #ifdef MOAB_HAVE_MPI
1165: /* filter all the non-owned and shared entities out of the list */
1166: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1167: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT); MBERRNM(merr);
1168: #endif
1170: /* get all the nodes via connectivity and the parent elements via adjacency information */
1171: merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false); MBERRNM(ierr);
1172: merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION); MBERRNM(ierr);
1173: }
1174: else {
1175: /* Let us query the hierarchy manager and get the results directly for this level */
1176: for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1177: moab::EntityHandle elemHandle = *iter;
1178: if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1179: dmmoab->bndyelems->insert(elemHandle);
1180: /* For this boundary element, query the vertices and add them to the list */
1181: std::vector<moab::EntityHandle> connect;
1182: merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect); MBERRNM(ierr);
1183: for (unsigned iv=0; iv < connect.size(); ++iv)
1184: if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv]))
1185: dmmoab->bndyvtx->insert(connect[iv]);
1186: /* Next, let us query the boundary faces and add them also to the list */
1187: std::vector<moab::EntityHandle> faces;
1188: merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim-1, faces); MBERRNM(ierr);
1189: for (unsigned ifa=0; ifa < faces.size(); ++ifa)
1190: if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa]))
1191: dmmoab->bndyfaces->insert(faces[ifa]);
1192: }
1193: }
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->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1197: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1198: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1199: #endif
1201: }
1202: PetscInfo3(NULL, "Found %D boundary vertices, %D boundary faces and %D boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size());
1204: /* Get the material sets and populate the data for all locally owned elements */
1205: {
1206: PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials);
1207: /* Get the count of entities of particular type from dmmoab->elocal
1208: -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1209: moab::Range msets;
1210: 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);
1211: if (msets.size() == 0) {
1212: PetscInfo(NULL, "No material sets found in the fileset.");
1213: }
1215: for (unsigned i=0; i < msets.size(); ++i) {
1216: moab::Range msetelems;
1217: merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);MB_CHK_ERR(merr);
1218: #ifdef MOAB_HAVE_MPI
1219: /* filter all the non-owned and shared entities out of the list */
1220: merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1221: #endif
1223: int partID;
1224: moab::EntityHandle mset=msets[i];
1225: merr = dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID);MB_CHK_ERR(merr);
1227: for (unsigned j=0; j < msetelems.size(); ++j)
1228: dmmoab->materials[dmmoab->elocal->index(msetelems[j])]=partID;
1229: }
1230: }
1232: return(0);
1233: }
1235: /*@C
1236: DMMoabCreateVertices - Creates and adds several vertices to the primary set represented by the DM.
1238: Collective
1240: Input Parameters:
1241: + dm - The DM object
1242: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1243: . conn - The connectivity of the element
1244: - nverts - The number of vertices that form the element
1246: Output Parameter:
1247: . overts - The list of vertices that were created (can be NULL)
1249: Level: beginner
1251: .seealso: DMMoabCreateSubmesh(), DMMoabCreateElement()
1252: @*/
1253: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal* coords, PetscInt nverts, moab::Range* overts)
1254: {
1255: moab::ErrorCode merr;
1256: DM_Moab *dmmoab;
1257: moab::Range verts;
1263: dmmoab = (DM_Moab*) dm->data;
1265: /* Insert new points */
1266: merr = dmmoab->mbiface->create_vertices(&coords[0], nverts, verts); MBERRNM(merr);
1267: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, verts); MBERRNM(merr);
1269: if (overts) *overts = verts;
1270: return(0);
1271: }
1273: /*@C
1274: DMMoabCreateElement - Adds an element of specified type to the primary set represented by the DM.
1276: Collective
1278: Input Parameters:
1279: + dm - The DM object
1280: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1281: . conn - The connectivity of the element
1282: - nverts - The number of vertices that form the element
1284: Output Parameter:
1285: . oelem - The handle to the element created and added to the DM object
1287: Level: beginner
1289: .seealso: DMMoabCreateSubmesh(), DMMoabCreateVertices()
1290: @*/
1291: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle* conn, PetscInt nverts, moab::EntityHandle* oelem)
1292: {
1293: moab::ErrorCode merr;
1294: DM_Moab *dmmoab;
1295: moab::EntityHandle elem;
1301: dmmoab = (DM_Moab*) dm->data;
1303: /* Insert new element */
1304: merr = dmmoab->mbiface->create_element(type, conn, nverts, elem); MBERRNM(merr);
1305: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1); MBERRNM(merr);
1307: if (oelem) *oelem = elem;
1308: return(0);
1309: }
1311: /*@C
1312: DMMoabCreateSubmesh - Creates a sub-DM object with a set that contains all vertices/elements of the parent
1313: in addition to providing support for dynamic mesh modifications. This is useful for AMR calculations to
1314: create a DM object on a refined level.
1316: Collective
1318: Input Parameters:
1319: . dm - The DM object
1321: Output Parameter:
1322: . newdm - The sub DM object with updated set information
1324: Level: advanced
1326: .seealso: DMCreate(), DMMoabCreateVertices(), DMMoabCreateElement()
1327: @*/
1328: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1329: {
1330: DM_Moab *dmmoab;
1331: DM_Moab *ndmmoab;
1332: moab::ErrorCode merr;
1333: PetscErrorCode ierr;
1338: dmmoab = (DM_Moab*) dm->data;
1340: /* Create the basic DMMoab object and keep the default parameters created by DM impls */
1341: DMMoabCreateMoab(((PetscObject)dm)->comm, dmmoab->mbiface, &dmmoab->ltog_tag, PETSC_NULL, newdm);
1343: /* get all the necessary handles from the private DM object */
1344: ndmmoab = (DM_Moab*) (*newdm)->data;
1346: /* set the sub-mesh's parent DM reference */
1347: ndmmoab->parent = &dm;
1349: /* create a file set to associate all entities in current mesh */
1350: merr = ndmmoab->mbiface->create_meshset(moab::MESHSET_SET, ndmmoab->fileset); MBERR("Creating file set failed", merr);
1352: /* create a meshset and then add old fileset as child */
1353: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal); MBERR("Adding child vertices to parent failed", merr);
1354: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal); MBERR("Adding child elements to parent failed", merr);
1356: /* preserve the field association between the parent and sub-mesh objects */
1357: DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames);
1358: return(0);
1359: }
1361: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1362: {
1363: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
1364: const char *name;
1365: MPI_Comm comm;
1366: PetscMPIInt size;
1367: PetscErrorCode ierr;
1370: PetscObjectGetComm((PetscObject)dm, &comm);
1371: MPI_Comm_size(comm, &size);
1372: PetscObjectGetName((PetscObject) dm, &name);
1373: PetscViewerASCIIPushTab(viewer);
1374: if (name) {PetscViewerASCIIPrintf(viewer, "%s in %D dimensions:\n", name, dmmoab->dim);}
1375: else {PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dmmoab->dim);}
1376: /* print details about the global mesh */
1377: {
1378: PetscViewerASCIIPushTab(viewer);
1379: PetscViewerASCIIPrintf(viewer, "Sizes: cells=%D, vertices=%D, blocks=%D\n", dmmoab->nele, dmmoab->n, dmmoab->bs);
1380: /* print boundary data */
1381: PetscViewerASCIIPrintf(viewer, "Boundary trace:\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size());
1382: {
1383: PetscViewerASCIIPushTab(viewer);
1384: PetscViewerASCIIPrintf(viewer, "cells=%D, faces=%D, vertices=%D\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size());
1385: PetscViewerASCIIPopTab(viewer);
1386: }
1387: /* print field data */
1388: PetscViewerASCIIPrintf(viewer, "Fields: %D components\n", dmmoab->numFields);
1389: {
1390: PetscViewerASCIIPushTab(viewer);
1391: for (int i = 0; i < dmmoab->numFields; ++i) {
1392: PetscViewerASCIIPrintf(viewer, "[%D] - %s\n", i, dmmoab->fieldNames[i]);
1393: }
1394: PetscViewerASCIIPopTab(viewer);
1395: }
1396: PetscViewerASCIIPopTab(viewer);
1397: }
1398: PetscViewerASCIIPopTab(viewer);
1399: PetscViewerFlush(viewer);
1400: return(0);
1401: }
1403: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1404: {
1405: return(0);
1406: }
1408: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1409: {
1410: return(0);
1411: }
1413: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1414: {
1415: PetscBool iascii, ishdf5, isvtk;
1421: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1422: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
1423: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
1424: if (iascii) {
1425: DMMoabView_Ascii(dm, viewer);
1426: } else if (ishdf5) {
1427: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1428: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
1429: DMMoabView_HDF5(dm, viewer);
1430: PetscViewerPopFormat(viewer);
1431: #else
1432: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1433: #endif
1434: }
1435: else if (isvtk) {
1436: DMMoabView_VTK(dm, viewer);
1437: }
1438: return(0);
1439: }
1441: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
1442: {
1444: dm->ops->view = DMView_Moab;
1445: dm->ops->load = NULL /* DMLoad_Moab */;
1446: dm->ops->setfromoptions = DMSetFromOptions_Moab;
1447: dm->ops->clone = DMClone_Moab;
1448: dm->ops->setup = DMSetUp_Moab;
1449: dm->ops->createlocalsection = NULL;
1450: dm->ops->createdefaultconstraints = NULL;
1451: dm->ops->createglobalvector = DMCreateGlobalVector_Moab;
1452: dm->ops->createlocalvector = DMCreateLocalVector_Moab;
1453: dm->ops->getlocaltoglobalmapping = NULL;
1454: dm->ops->createfieldis = NULL;
1455: dm->ops->createcoordinatedm = NULL /* DMCreateCoordinateDM_Moab */;
1456: dm->ops->getcoloring = NULL;
1457: dm->ops->creatematrix = DMCreateMatrix_Moab;
1458: dm->ops->createinterpolation = DMCreateInterpolation_Moab;
1459: dm->ops->createinjection = NULL /* DMCreateInjection_Moab */;
1460: dm->ops->refine = DMRefine_Moab;
1461: dm->ops->coarsen = DMCoarsen_Moab;
1462: dm->ops->refinehierarchy = DMRefineHierarchy_Moab;
1463: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Moab;
1464: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Moab;
1465: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Moab;
1466: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Moab;
1467: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Moab;
1468: dm->ops->destroy = DMDestroy_Moab;
1469: dm->ops->createsubdm = NULL /* DMCreateSubDM_Moab */;
1470: dm->ops->getdimpoints = NULL /* DMGetDimPoints_Moab */;
1471: dm->ops->locatepoints = NULL /* DMLocatePoints_Moab */;
1472: return(0);
1473: }
1475: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1476: {
1477: PetscErrorCode ierr;
1480: /* get all the necessary handles from the private DM object */
1481: (*newdm)->data = (DM_Moab*) dm->data;
1482: ((DM_Moab*)dm->data)->refct++;
1484: PetscObjectChangeTypeName((PetscObject) *newdm, DMMOAB);
1485: DMInitialize_Moab(*newdm);
1486: return(0);
1487: }
1489: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1490: {
1495: PetscNewLog(dm, (DM_Moab**)&dm->data);
1497: ((DM_Moab*)dm->data)->bs = 1;
1498: ((DM_Moab*)dm->data)->numFields = 1;
1499: ((DM_Moab*)dm->data)->n = 0;
1500: ((DM_Moab*)dm->data)->nloc = 0;
1501: ((DM_Moab*)dm->data)->nghost = 0;
1502: ((DM_Moab*)dm->data)->nele = 0;
1503: ((DM_Moab*)dm->data)->neleloc = 0;
1504: ((DM_Moab*)dm->data)->neleghost = 0;
1505: ((DM_Moab*)dm->data)->ltog_map = NULL;
1506: ((DM_Moab*)dm->data)->ltog_sendrecv = NULL;
1508: ((DM_Moab*)dm->data)->refct = 1;
1509: ((DM_Moab*)dm->data)->parent = NULL;
1510: ((DM_Moab*)dm->data)->vlocal = new moab::Range();
1511: ((DM_Moab*)dm->data)->vowned = new moab::Range();
1512: ((DM_Moab*)dm->data)->vghost = new moab::Range();
1513: ((DM_Moab*)dm->data)->elocal = new moab::Range();
1514: ((DM_Moab*)dm->data)->eghost = new moab::Range();
1516: DMInitialize_Moab(dm);
1517: return(0);
1518: }