Actual source code: dmmoab.cxx
1: #include <petsc/private/dmmbimpl.h>
3: #include <petscdmmoab.h>
4: #include <MBTagConventions.hpp>
5: #include <moab/NestedRefine.hpp>
6: #include <moab/Skinner.hpp>
8: /*MC
9: DMMOAB = "moab" - A DM object that encapsulates an unstructured mesh described by the MOAB mesh database.
10: Direct access to the MOAB Interface and other mesh manipulation related objects are available
11: through public API. Ability to create global and local representation of Vecs containing all
12: unknowns in the interior and shared boundary via a transparent tag-data wrapper is provided
13: along with utility functions to traverse the mesh and assemble a discrete system via
14: field-based/blocked Vec(Get/Set) methods. Input from and output to different formats are
15: available.
17: Reference: https://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html
19: Level: intermediate
21: .seealso: DMType, DMMoabCreate(), DMCreate(), DMSetType(), DMMoabCreateMoab()
22: M*/
24: /* External function declarations here */
25: PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
26: PETSC_EXTERN PetscErrorCode DMCreateDefaultConstraints_Moab(DM dm);
27: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
28: PETSC_EXTERN PetscErrorCode DMCreateCoordinateDM_Moab(DM dm, DM *cdm);
29: PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM *dmRefined);
30: PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM *dmCoarsened);
31: PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmRefined[]);
32: PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmCoarsened[]);
33: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm);
34: PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM, Vec *);
35: PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM, Vec *);
36: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
37: PETSC_EXTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM, Vec, InsertMode, Vec);
38: PETSC_EXTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM, Vec, InsertMode, Vec);
39: PETSC_EXTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM, Vec, InsertMode, Vec);
40: PETSC_EXTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM, Vec, InsertMode, Vec);
42: /* Un-implemented routines */
43: /*
44: PETSC_EXTERN PetscErrorCode DMCreatelocalsection_Moab(DM dm);
45: PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dmCoarse, DM dmFine, Mat *mat);
46: PETSC_EXTERN PetscErrorCode DMLoad_Moab(DM dm, PetscViewer viewer);
47: PETSC_EXTERN PetscErrorCode DMGetDimPoints_Moab(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd);
48: PETSC_EXTERN PetscErrorCode DMCreateSubDM_Moab(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
49: PETSC_EXTERN PetscErrorCode DMLocatePoints_Moab(DM dm, Vec v, IS *cellIS);
50: */
52: /*@C
53: DMMoabCreate - Creates a DMMoab object, which encapsulates a moab instance
55: Collective
57: Input Parameter:
58: . comm - The communicator for the DMMoab object
60: Output Parameter:
61: . dmb - The DMMoab object
63: Level: beginner
65: @*/
66: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
67: {
69: DMCreate(comm, dmb);
70: DMSetType(*dmb, DMMOAB);
71: return 0;
72: }
74: /*@C
75: DMMoabCreateMoab - Creates a DMMoab object, optionally from an instance and other data
77: Collective
79: Input Parameters:
80: + comm - The communicator for the DMMoab object
81: . mbiface - (ptr to) the MOAB Instance; if passed in NULL, MOAB instance is created inside PETSc, and destroyed
82: along with the DMMoab
83: . ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
84: - range - If non-NULL, contains range of entities to which DOFs will be assigned
86: Output Parameter:
87: . dmb - The DMMoab object
89: Level: intermediate
91: @*/
92: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
93: {
94: moab::ErrorCode merr;
95: DM dmmb;
96: DM_Moab *dmmoab;
100: DMMoabCreate(comm, &dmmb);
101: dmmoab = (DM_Moab*)(dmmb)->data;
103: if (!mbiface) {
104: dmmoab->mbiface = new moab::Core();
105: dmmoab->icreatedinstance = PETSC_TRUE;
106: }
107: else {
108: dmmoab->mbiface = mbiface;
109: dmmoab->icreatedinstance = PETSC_FALSE;
110: }
112: /* by default the fileset = root set. This set stores the hierarchy of entities belonging to current DM */
113: dmmoab->fileset = 0;
114: dmmoab->hlevel = 0;
115: dmmoab->nghostrings = 0;
117: #ifdef MOAB_HAVE_MPI
118: moab::EntityHandle partnset;
120: /* Create root sets for each mesh. Then pass these
121: to the load_file functions to be populated. */
122: merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, partnset); MBERR("Creating partition set failed", merr);
124: /* Create the parallel communicator object with the partition handle associated with MOAB */
125: dmmoab->pcomm = moab::ParallelComm::get_pcomm(dmmoab->mbiface, partnset, &comm);
126: #endif
128: /* do the remaining initializations for DMMoab */
129: dmmoab->bs = 1;
130: dmmoab->numFields = 1;
131: PetscMalloc(dmmoab->numFields * sizeof(char*), &dmmoab->fieldNames);
132: PetscStrallocpy("DEFAULT", (char**) &dmmoab->fieldNames[0]);
133: dmmoab->rw_dbglevel = 0;
134: dmmoab->partition_by_rank = PETSC_FALSE;
135: dmmoab->extra_read_options[0] = '\0';
136: dmmoab->extra_write_options[0] = '\0';
137: dmmoab->read_mode = READ_PART;
138: dmmoab->write_mode = WRITE_PART;
140: /* set global ID tag handle */
141: if (ltog_tag && *ltog_tag) {
142: DMMoabSetLocalToGlobalTag(dmmb, *ltog_tag);
143: }
144: else {
145: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag); MBERRNM(merr);
146: if (ltog_tag) *ltog_tag = dmmoab->ltog_tag;
147: }
149: merr = dmmoab->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dmmoab->material_tag); MBERRNM(merr);
151: /* set the local range of entities (vertices) of interest */
152: if (range) {
153: DMMoabSetLocalVertices(dmmb, range);
154: }
155: *dmb = dmmb;
156: return 0;
157: }
159: #ifdef MOAB_HAVE_MPI
161: /*@C
162: DMMoabGetParallelComm - Get the ParallelComm used with this DMMoab
164: Collective
166: Input Parameter:
167: . dm - The DMMoab object being set
169: Output Parameter:
170: . pcomm - The ParallelComm for the DMMoab
172: Level: beginner
174: @*/
175: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
176: {
178: *pcomm = ((DM_Moab*)(dm)->data)->pcomm;
179: return 0;
180: }
182: #endif /* MOAB_HAVE_MPI */
184: /*@C
185: DMMoabSetInterface - Set the MOAB instance used with this DMMoab
187: Collective
189: Input Parameters:
190: + dm - The DMMoab object being set
191: - mbiface - The MOAB instance being set on this DMMoab
193: Level: beginner
195: @*/
196: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
197: {
198: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
202: #ifdef MOAB_HAVE_MPI
203: dmmoab->pcomm = NULL;
204: #endif
205: dmmoab->mbiface = mbiface;
206: dmmoab->icreatedinstance = PETSC_FALSE;
207: return 0;
208: }
210: /*@C
211: DMMoabGetInterface - Get the MOAB instance used with this DMMoab
213: Collective
215: Input Parameter:
216: . dm - The DMMoab object being set
218: Output Parameter:
219: . mbiface - The MOAB instance set on this DMMoab
221: Level: beginner
223: @*/
224: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
225: {
226: static PetscBool cite = PETSC_FALSE;
229: 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);
230: *mbiface = ((DM_Moab*)dm->data)->mbiface;
231: return 0;
232: }
234: /*@C
235: DMMoabSetLocalVertices - Set the entities having DOFs on this DMMoab
237: Collective
239: Input Parameters:
240: + dm - The DMMoab object being set
241: - range - The entities treated by this DMMoab
243: Level: beginner
245: @*/
246: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
247: {
248: moab::Range tmpvtxs;
249: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
252: dmmoab->vlocal->clear();
253: dmmoab->vowned->clear();
255: dmmoab->vlocal->insert(range->begin(), range->end());
257: #ifdef MOAB_HAVE_MPI
258: moab::ErrorCode merr;
259: /* filter based on parallel status */
260: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned); MBERRNM(merr);
262: /* filter all the non-owned and shared entities out of the list */
263: tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
264: merr = dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost); MBERRNM(merr);
265: tmpvtxs = moab::subtract(tmpvtxs, *dmmoab->vghost);
266: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
267: #else
268: *dmmoab->vowned = *dmmoab->vlocal;
269: #endif
271: /* compute and cache the sizes of local and ghosted entities */
272: dmmoab->nloc = dmmoab->vowned->size();
273: dmmoab->nghost = dmmoab->vghost->size();
274: #ifdef MOAB_HAVE_MPI
275: MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
276: #else
277: dmmoab->n = dmmoab->nloc;
278: #endif
279: return 0;
280: }
282: /*@C
283: DMMoabGetAllVertices - Get the entities having DOFs on this DMMoab
285: Collective
287: Input Parameter:
288: . dm - The DMMoab object being set
290: Output Parameter:
291: . owned - The local vertex entities in this DMMoab = (owned+ghosted)
293: Level: beginner
295: @*/
296: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
297: {
299: if (local) *local = *((DM_Moab*)dm->data)->vlocal;
300: return 0;
301: }
303: /*@C
304: DMMoabGetLocalVertices - Get the entities having DOFs on this DMMoab
306: Collective
308: Input Parameter:
309: . dm - The DMMoab object being set
311: Output Parameters:
312: + owned - The owned vertex entities in this DMMoab
313: - ghost - The ghosted entities (non-owned) stored locally in this partition
315: Level: beginner
317: @*/
318: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
319: {
321: if (owned) *owned = ((DM_Moab*)dm->data)->vowned;
322: if (ghost) *ghost = ((DM_Moab*)dm->data)->vghost;
323: return 0;
324: }
326: /*@C
327: DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned
329: Collective
331: Input Parameter:
332: . dm - The DMMoab object being set
334: Output Parameter:
335: . range - The entities owned locally
337: Level: beginner
339: @*/
340: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
341: {
343: if (range) *range = ((DM_Moab*)dm->data)->elocal;
344: return 0;
345: }
347: /*@C
348: DMMoabSetLocalElements - Set the entities having DOFs on this DMMoab
350: Collective
352: Input Parameters:
353: + dm - The DMMoab object being set
354: - range - The entities treated by this DMMoab
356: Level: beginner
358: @*/
359: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
360: {
361: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
364: dmmoab->elocal->clear();
365: dmmoab->eghost->clear();
366: dmmoab->elocal->insert(range->begin(), range->end());
367: #ifdef MOAB_HAVE_MPI
368: moab::ErrorCode merr;
369: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
370: *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
371: #endif
372: dmmoab->neleloc = dmmoab->elocal->size();
373: dmmoab->neleghost = dmmoab->eghost->size();
374: #ifdef MOAB_HAVE_MPI
375: MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
376: PetscInfo(dm, "Created %D local and %D global elements.\n", dmmoab->neleloc, dmmoab->nele);
377: #else
378: dmmoab->nele = dmmoab->neleloc;
379: #endif
380: return 0;
381: }
383: /*@C
384: DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
386: Collective
388: Input Parameters:
389: + dm - The DMMoab object being set
390: - ltogtag - The MOAB tag used for local to global ids
392: Level: beginner
394: @*/
395: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
396: {
398: ((DM_Moab*)dm->data)->ltog_tag = ltogtag;
399: return 0;
400: }
402: /*@C
403: DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
405: Collective
407: Input Parameter:
408: . dm - The DMMoab object being set
410: Output Parameter:
411: . ltogtag - The MOAB tag used for local to global ids
413: Level: beginner
415: @*/
416: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
417: {
419: *ltog_tag = ((DM_Moab*)dm->data)->ltog_tag;
420: return 0;
421: }
423: /*@C
424: DMMoabSetBlockSize - Set the block size used with this DMMoab
426: Collective
428: Input Parameters:
429: + dm - The DMMoab object being set
430: - bs - The block size used with this DMMoab
432: Level: beginner
434: @*/
435: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
436: {
438: ((DM_Moab*)dm->data)->bs = bs;
439: return 0;
440: }
442: /*@C
443: DMMoabGetBlockSize - Get the block size used with this DMMoab
445: Collective
447: Input Parameter:
448: . dm - The DMMoab object being set
450: Output Parameter:
451: . bs - The block size used with this DMMoab
453: Level: beginner
455: @*/
456: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
457: {
459: *bs = ((DM_Moab*)dm->data)->bs;
460: return 0;
461: }
463: /*@C
464: DMMoabGetSize - Get the global vertex size used with this DMMoab
466: Collective on dm
468: Input Parameter:
469: . dm - The DMMoab object being set
471: Output Parameters:
472: + neg - The number of global elements in the DMMoab instance
473: - nvg - The number of global vertices in the DMMoab instance
475: Level: beginner
477: @*/
478: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
479: {
481: if (neg) *neg = ((DM_Moab*)dm->data)->nele;
482: if (nvg) *nvg = ((DM_Moab*)dm->data)->n;
483: return 0;
484: }
486: /*@C
487: DMMoabGetLocalSize - Get the local and ghosted vertex size used with this DMMoab
489: Collective on dm
491: Input Parameter:
492: . dm - The DMMoab object being set
494: Output Parameters:
495: + nel - The number of owned elements in this processor
496: . neg - The number of ghosted elements in this processor
497: . nvl - The number of owned vertices in this processor
498: - nvg - The number of ghosted vertices in this processor
500: Level: beginner
502: @*/
503: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
504: {
506: if (nel) *nel = ((DM_Moab*)dm->data)->neleloc;
507: if (neg) *neg = ((DM_Moab*)dm->data)->neleghost;
508: if (nvl) *nvl = ((DM_Moab*)dm->data)->nloc;
509: if (nvg) *nvg = ((DM_Moab*)dm->data)->nghost;
510: return 0;
511: }
513: /*@C
514: DMMoabGetOffset - Get the local offset for the global vector
516: Collective
518: Input Parameter:
519: . dm - The DMMoab object being set
521: Output Parameter:
522: . offset - The local offset for the global vector
524: Level: beginner
526: @*/
527: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
528: {
530: *offset = ((DM_Moab*)dm->data)->vstart;
531: return 0;
532: }
534: /*@C
535: DMMoabGetDimension - Get the dimension of the DM Mesh
537: Collective
539: Input Parameter:
540: . dm - The DMMoab object
542: Output Parameter:
543: . dim - The dimension of DM
545: Level: beginner
547: @*/
548: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
549: {
551: *dim = ((DM_Moab*)dm->data)->dim;
552: return 0;
553: }
555: /*@C
556: DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
557: generated through uniform refinement.
559: Collective on dm
561: Input Parameter:
562: . dm - The DMMoab object being set
564: Output Parameter:
565: . nvg - The current mesh hierarchy level
567: Level: beginner
569: @*/
570: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
571: {
573: if (nlevel) *nlevel = ((DM_Moab*)dm->data)->hlevel;
574: return 0;
575: }
577: /*@C
578: DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh
580: Collective
582: Input Parameters:
583: + dm - The DMMoab object
584: - ehandle - The element entity handle
586: Output Parameter:
587: . mat - The material ID for the current entity
589: Level: beginner
591: @*/
592: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
593: {
594: DM_Moab *dmmoab;
597: if (*mat) {
598: dmmoab = (DM_Moab*)(dm)->data;
599: *mat = dmmoab->materials[dmmoab->elocal->index(ehandle)];
600: }
601: return 0;
602: }
604: /*@C
605: DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities
607: Collective
609: Input Parameters:
610: + dm - The DMMoab object
611: . nconn - Number of entities whose coordinates are needed
612: - conn - The vertex entity handles
614: Output Parameter:
615: . vpos - The coordinates of the requested vertex entities
617: Level: beginner
619: .seealso: DMMoabGetVertexConnectivity()
620: @*/
621: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
622: {
623: DM_Moab *dmmoab;
624: moab::ErrorCode merr;
629: dmmoab = (DM_Moab*)(dm)->data;
631: /* Get connectivity information in MOAB canonical ordering */
632: if (dmmoab->hlevel) {
633: merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle*>(conn), nconn, dmmoab->hlevel, vpos);MBERRNM(merr);
634: }
635: else {
636: merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);MBERRNM(merr);
637: }
638: return 0;
639: }
641: /*@C
642: DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity
644: Collective
646: Input Parameters:
647: + dm - The DMMoab object
648: - vhandle - Vertex entity handle
650: Output Parameters:
651: + nconn - Number of entities whose coordinates are needed
652: - conn - The vertex entity handles
654: Level: beginner
656: .seealso: DMMoabGetVertexCoordinates(), DMMoabRestoreVertexConnectivity()
657: @*/
658: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt* nconn, moab::EntityHandle **conn)
659: {
660: DM_Moab *dmmoab;
661: std::vector<moab::EntityHandle> adj_entities, connect;
662: moab::ErrorCode merr;
666: dmmoab = (DM_Moab*)(dm)->data;
668: /* Get connectivity information in MOAB canonical ordering */
669: merr = dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION); MBERRNM(merr);
670: merr = dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect); MBERRNM(merr);
672: if (conn) {
673: PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn);
674: PetscArraycpy(*conn, &connect[0], connect.size());
675: }
676: if (nconn) *nconn = connect.size();
677: return 0;
678: }
680: /*@C
681: DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity
683: Collective
685: Input Parameters:
686: + dm - The DMMoab object
687: . vhandle - Vertex entity handle
688: . nconn - Number of entities whose coordinates are needed
689: - conn - The vertex entity handles
691: Level: beginner
693: .seealso: DMMoabGetVertexCoordinates(), DMMoabGetVertexConnectivity()
694: @*/
695: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt* nconn, moab::EntityHandle **conn)
696: {
700: if (conn) {
701: PetscFree(*conn);
702: }
703: if (nconn) *nconn = 0;
704: return 0;
705: }
707: /*@C
708: DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity
710: Collective
712: Input Parameters:
713: + dm - The DMMoab object
714: - ehandle - Vertex entity handle
716: Output Parameters:
717: + nconn - Number of entities whose coordinates are needed
718: - conn - The vertex entity handles
720: Level: beginner
722: .seealso: DMMoabGetVertexCoordinates(), DMMoabGetVertexConnectivity(), DMMoabRestoreVertexConnectivity()
723: @*/
724: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt* nconn, const moab::EntityHandle **conn)
725: {
726: DM_Moab *dmmoab;
727: const moab::EntityHandle *connect;
728: std::vector<moab::EntityHandle> vconn;
729: moab::ErrorCode merr;
730: PetscInt nnodes;
734: dmmoab = (DM_Moab*)(dm)->data;
736: /* Get connectivity information in MOAB canonical ordering */
737: merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes); MBERRNM(merr);
738: if (conn) *conn = connect;
739: if (nconn) *nconn = nnodes;
740: return 0;
741: }
743: /*@C
744: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
746: Collective
748: Input Parameters:
749: + dm - The DMMoab object
750: - ent - Entity handle
752: Output Parameter:
753: . ent_on_boundary - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
755: Level: beginner
757: .seealso: DMMoabCheckBoundaryVertices()
758: @*/
759: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool* ent_on_boundary)
760: {
761: moab::EntityType etype;
762: DM_Moab *dmmoab;
763: PetscInt edim;
767: dmmoab = (DM_Moab*)(dm)->data;
769: /* get the entity type and handle accordingly */
770: etype = dmmoab->mbiface->type_from_handle(ent);
773: /* get the entity dimension */
774: edim = dmmoab->mbiface->dimension_from_handle(ent);
776: *ent_on_boundary = PETSC_FALSE;
777: if (etype == moab::MBVERTEX && edim == 0) {
778: *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
779: }
780: else {
781: if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
782: if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
783: }
784: else { /* next check the lower-dimensional faces */
785: if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
786: }
787: }
788: return 0;
789: }
791: /*@C
792: DMMoabCheckBoundaryVertices - Check whether a given entity is on the boundary (vertex, edge, face, element)
794: Input Parameters:
795: + dm - The DMMoab object
796: . nconn - Number of handles
797: - cnt - Array of entity handles
799: Output Parameter:
800: . isbdvtx - Array of boundary markers - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
802: Level: beginner
804: .seealso: DMMoabIsEntityOnBoundary()
805: @*/
806: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool* isbdvtx)
807: {
808: DM_Moab *dmmoab;
809: PetscInt i;
814: dmmoab = (DM_Moab*)(dm)->data;
816: for (i = 0; i < nconn; ++i) {
817: isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
818: }
819: return 0;
820: }
822: /*@C
823: DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary
825: Input Parameter:
826: . dm - The DMMoab object
828: Output Parameters:
829: + bdvtx - Boundary vertices
830: . bdelems - Boundary elements
831: - bdfaces - Boundary faces
833: Level: beginner
835: .seealso: DMMoabCheckBoundaryVertices(), DMMoabIsEntityOnBoundary()
836: @*/
837: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range** bdelems, const moab::Range** bdfaces)
838: {
839: DM_Moab *dmmoab;
842: dmmoab = (DM_Moab*)(dm)->data;
844: if (bdvtx) *bdvtx = dmmoab->bndyvtx;
845: if (bdfaces) *bdfaces = dmmoab->bndyfaces;
846: if (bdelems) *bdfaces = dmmoab->bndyelems;
847: return 0;
848: }
850: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
851: {
852: PetscInt i;
853: moab::ErrorCode merr;
854: DM_Moab *dmmoab = (DM_Moab*)dm->data;
858: dmmoab->refct--;
859: if (!dmmoab->refct) {
860: delete dmmoab->vlocal;
861: delete dmmoab->vowned;
862: delete dmmoab->vghost;
863: delete dmmoab->elocal;
864: delete dmmoab->eghost;
865: delete dmmoab->bndyvtx;
866: delete dmmoab->bndyfaces;
867: delete dmmoab->bndyelems;
869: PetscFree(dmmoab->gsindices);
870: PetscFree2(dmmoab->gidmap, dmmoab->lidmap);
871: PetscFree(dmmoab->dfill);
872: PetscFree(dmmoab->ofill);
873: PetscFree(dmmoab->materials);
874: if (dmmoab->fieldNames) {
875: for (i = 0; i < dmmoab->numFields; i++) {
876: PetscFree(dmmoab->fieldNames[i]);
877: }
878: PetscFree(dmmoab->fieldNames);
879: }
881: if (dmmoab->nhlevels) {
882: PetscFree(dmmoab->hsets);
883: dmmoab->nhlevels = 0;
884: if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
885: dmmoab->hierarchy = NULL;
886: }
888: if (dmmoab->icreatedinstance) {
889: delete dmmoab->pcomm;
890: merr = dmmoab->mbiface->delete_mesh(); MBERRNM(merr);
891: delete dmmoab->mbiface;
892: }
893: dmmoab->mbiface = NULL;
894: #ifdef MOAB_HAVE_MPI
895: dmmoab->pcomm = NULL;
896: #endif
897: VecScatterDestroy(&dmmoab->ltog_sendrecv);
898: ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map);
899: PetscFree(dm->data);
900: }
901: return 0;
902: }
904: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(PetscOptionItems *PetscOptionsObject, DM dm)
905: {
906: DM_Moab *dmmoab = (DM_Moab*)dm->data;
909: PetscOptionsHead(PetscOptionsObject, "DMMoab Options");
910: PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL,0);
911: PetscOptionsBool("-dm_moab_partiton_by_rank", "Use partition by rank when reading MOAB meshes from file", "DMView", dmmoab->partition_by_rank, &dmmoab->partition_by_rank, NULL);
912: /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
913: PetscOptionsString("-dm_moab_read_opts", "Extra options to enable MOAB reader to load DM from file", "DMView", dmmoab->extra_read_options, dmmoab->extra_read_options, sizeof(dmmoab->extra_read_options), NULL);
914: PetscOptionsString("-dm_moab_write_opts", "Extra options to enable MOAB writer to serialize DM to file", "DMView", dmmoab->extra_write_options, dmmoab->extra_write_options, sizeof(dmmoab->extra_write_options), NULL);
915: PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum*)&dmmoab->read_mode, NULL);
916: PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum*)&dmmoab->write_mode, NULL);
917: return 0;
918: }
920: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
921: {
922: moab::ErrorCode merr;
923: Vec local, global;
924: IS from, to;
925: moab::Range::iterator iter;
926: PetscInt i, j, f, bs, vent, totsize, *lgmap;
927: DM_Moab *dmmoab = (DM_Moab*)dm->data;
928: moab::Range adjs;
931: /* Get the local and shared vertices and cache it */
933: #ifdef MOAB_HAVE_MPI
935: #endif
937: /* Get the entities recursively in the current part of the mesh, if user did not set the local vertices explicitly */
938: if (dmmoab->vlocal->empty())
939: {
940: //merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,*dmmoab->vlocal,true);MBERRNM(merr);
941: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, 0, *dmmoab->vlocal, false); MBERRNM(merr);
943: #ifdef MOAB_HAVE_MPI
944: /* filter based on parallel status */
945: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned); MBERRNM(merr);
947: /* filter all the non-owned and shared entities out of the list */
948: // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
949: adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
950: merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost); MBERRNM(merr);
951: adjs = moab::subtract(adjs, *dmmoab->vghost);
952: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
953: #else
954: *dmmoab->vowned = *dmmoab->vlocal;
955: #endif
957: /* compute and cache the sizes of local and ghosted entities */
958: dmmoab->nloc = dmmoab->vowned->size();
959: dmmoab->nghost = dmmoab->vghost->size();
961: #ifdef MOAB_HAVE_MPI
962: MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
963: PetscInfo(NULL, "Filset ID: %u, Vertices: local - %D, owned - %D, ghosted - %D.\n", dmmoab->fileset, dmmoab->vlocal->size(), dmmoab->nloc, dmmoab->nghost);
964: #else
965: dmmoab->n = dmmoab->nloc;
966: #endif
967: }
969: {
970: /* get the information about the local elements in the mesh */
971: dmmoab->eghost->clear();
973: /* first decipher the leading dimension */
974: for (i = 3; i > 0; i--) {
975: dmmoab->elocal->clear();
976: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false); MBERRNM(merr);
978: /* store the current mesh dimension */
979: if (dmmoab->elocal->size()) {
980: dmmoab->dim = i;
981: break;
982: }
983: }
985: DMSetDimension(dm, dmmoab->dim);
987: #ifdef MOAB_HAVE_MPI
988: /* filter the ghosted and owned element list */
989: *dmmoab->eghost = *dmmoab->elocal;
990: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
991: *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
992: #endif
994: dmmoab->neleloc = dmmoab->elocal->size();
995: dmmoab->neleghost = dmmoab->eghost->size();
997: #ifdef MOAB_HAVE_MPI
998: MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
999: PetscInfo(NULL, "%d-dim elements: owned - %D, ghosted - %D.\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost);
1000: #else
1001: dmmoab->nele = dmmoab->neleloc;
1002: #endif
1003: }
1005: bs = dmmoab->bs;
1006: if (!dmmoab->ltog_tag) {
1007: /* Get the global ID tag. The global ID tag is applied to each
1008: vertex. It acts as an global identifier which MOAB uses to
1009: assemble the individual pieces of the mesh */
1010: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag); MBERRNM(merr);
1011: }
1013: totsize = dmmoab->vlocal->size();
1015: PetscCalloc1(totsize, &dmmoab->gsindices);
1016: {
1017: /* first get the local indices */
1018: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]); MBERRNM(merr);
1019: if (dmmoab->nghost) { /* next get the ghosted indices */
1020: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]); MBERRNM(merr);
1021: }
1023: /* find out the local and global minima of GLOBAL_ID */
1024: dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1025: for (i = 0; i < totsize; ++i) {
1026: if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1027: if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1028: }
1030: MPIU_Allreduce(&dmmoab->lminmax[0], &dmmoab->gminmax[0], 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm);
1031: MPIU_Allreduce(&dmmoab->lminmax[1], &dmmoab->gminmax[1], 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm);
1033: /* set the GID map */
1034: for (i = 0; i < totsize; ++i) {
1035: dmmoab->gsindices[i] -= dmmoab->gminmax[0]; /* zero based index needed for IS */
1037: }
1038: dmmoab->lminmax[0] -= dmmoab->gminmax[0];
1039: dmmoab->lminmax[1] -= dmmoab->gminmax[0];
1041: PetscInfo(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]);
1042: }
1045: {
1046: dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1047: dmmoab->seqend = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1048: PetscInfo(NULL, "SEQUENCE: Local [min, max] - [%D, %D]\n", dmmoab->seqstart, dmmoab->seqend);
1050: PetscMalloc2(dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->gidmap, dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->lidmap);
1051: PetscMalloc1(totsize * dmmoab->numFields, &lgmap);
1053: i = j = 0;
1054: /* set the owned vertex data first */
1055: for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1056: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1057: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1058: dmmoab->lidmap[vent] = i;
1059: for (f = 0; f < dmmoab->numFields; f++, j++) {
1060: lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1061: }
1062: }
1063: /* next arrange all the ghosted data information */
1064: for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1065: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1066: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1067: dmmoab->lidmap[vent] = i;
1068: for (f = 0; f < dmmoab->numFields; f++, j++) {
1069: lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1070: }
1071: }
1073: /* We need to create the Global to Local Vector Scatter Contexts
1074: 1) First create a local and global vector
1075: 2) Create a local and global IS
1076: 3) Create VecScatter and LtoGMapping objects
1077: 4) Cleanup the IS and Vec objects
1078: */
1079: DMCreateGlobalVector(dm, &global);
1080: DMCreateLocalVector(dm, &local);
1082: VecGetOwnershipRange(global, &dmmoab->vstart, &dmmoab->vend);
1084: /* global to local must retrieve ghost points */
1085: ISCreateStride(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, dmmoab->vstart, 1, &from);
1086: ISSetBlockSize(from, bs);
1088: ISCreateGeneral(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, &lgmap[0], PETSC_COPY_VALUES, &to);
1089: ISSetBlockSize(to, bs);
1091: if (!dmmoab->ltog_map) {
1092: /* create to the local to global mapping for vectors in order to use VecSetValuesLocal */
1093: ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm, dmmoab->bs, totsize * dmmoab->numFields, lgmap,PETSC_COPY_VALUES, &dmmoab->ltog_map);
1094: }
1096: /* now create the scatter object from local to global vector */
1097: VecScatterCreate(local, from, global, to, &dmmoab->ltog_sendrecv);
1099: /* clean up IS, Vec */
1100: PetscFree(lgmap);
1101: ISDestroy(&from);
1102: ISDestroy(&to);
1103: VecDestroy(&local);
1104: VecDestroy(&global);
1105: }
1107: dmmoab->bndyvtx = new moab::Range();
1108: dmmoab->bndyfaces = new moab::Range();
1109: dmmoab->bndyelems = new moab::Range();
1110: /* skin the boundary and store nodes */
1111: if (!dmmoab->hlevel) {
1112: /* get the skin vertices of boundary faces for the current partition and then filter
1113: the local, boundary faces, vertices and elements alone via PSTATUS flags;
1114: this should not give us any ghosted boundary, but if user needs such a functionality
1115: it would be easy to add it based on the find_skin query below */
1116: moab::Skinner skinner(dmmoab->mbiface);
1118: /* get the entities on the skin - only the faces */
1119: 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
1121: #ifdef MOAB_HAVE_MPI
1122: /* filter all the non-owned and shared entities out of the list */
1123: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);MBERRNM(merr);
1124: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT);MBERRNM(merr);
1125: #endif
1127: /* get all the nodes via connectivity and the parent elements via adjacency information */
1128: merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false);MBERRNM(merr);
1129: merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION);MBERRNM(merr);
1130: }
1131: else {
1132: /* Let us query the hierarchy manager and get the results directly for this level */
1133: for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1134: moab::EntityHandle elemHandle = *iter;
1135: if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1136: dmmoab->bndyelems->insert(elemHandle);
1137: /* For this boundary element, query the vertices and add them to the list */
1138: std::vector<moab::EntityHandle> connect;
1139: merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect);MBERRNM(merr);
1140: for (unsigned iv=0; iv < connect.size(); ++iv)
1141: if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv]))
1142: dmmoab->bndyvtx->insert(connect[iv]);
1143: /* Next, let us query the boundary faces and add them also to the list */
1144: std::vector<moab::EntityHandle> faces;
1145: merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim-1, faces);MBERRNM(merr);
1146: for (unsigned ifa=0; ifa < faces.size(); ++ifa)
1147: if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa]))
1148: dmmoab->bndyfaces->insert(faces[ifa]);
1149: }
1150: }
1151: #ifdef MOAB_HAVE_MPI
1152: /* filter all the non-owned and shared entities out of the list */
1153: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1154: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1155: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1156: #endif
1158: }
1159: PetscInfo(NULL, "Found %D boundary vertices, %D boundary faces and %D boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size());
1161: /* Get the material sets and populate the data for all locally owned elements */
1162: {
1163: PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials);
1164: /* Get the count of entities of particular type from dmmoab->elocal
1165: -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1166: moab::Range msets;
1167: 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);
1168: if (msets.size() == 0) {
1169: PetscInfo(NULL, "No material sets found in the fileset.");
1170: }
1172: for (unsigned i=0; i < msets.size(); ++i) {
1173: moab::Range msetelems;
1174: merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);MB_CHK_ERR(merr);
1175: #ifdef MOAB_HAVE_MPI
1176: /* filter all the non-owned and shared entities out of the list */
1177: merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT); MBERRNM(merr);
1178: #endif
1180: int partID;
1181: moab::EntityHandle mset=msets[i];
1182: merr = dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID);MB_CHK_ERR(merr);
1184: for (unsigned j=0; j < msetelems.size(); ++j)
1185: dmmoab->materials[dmmoab->elocal->index(msetelems[j])]=partID;
1186: }
1187: }
1189: return 0;
1190: }
1192: /*@C
1193: DMMoabCreateVertices - Creates and adds several vertices to the primary set represented by the DM.
1195: Collective
1197: Input Parameters:
1198: + dm - The DM object
1199: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1200: . conn - The connectivity of the element
1201: - nverts - The number of vertices that form the element
1203: Output Parameter:
1204: . overts - The list of vertices that were created (can be NULL)
1206: Level: beginner
1208: .seealso: DMMoabCreateSubmesh(), DMMoabCreateElement()
1209: @*/
1210: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal* coords, PetscInt nverts, moab::Range* overts)
1211: {
1212: moab::ErrorCode merr;
1213: DM_Moab *dmmoab;
1214: moab::Range verts;
1219: dmmoab = (DM_Moab*) dm->data;
1221: /* Insert new points */
1222: merr = dmmoab->mbiface->create_vertices(&coords[0], nverts, verts); MBERRNM(merr);
1223: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, verts); MBERRNM(merr);
1225: if (overts) *overts = verts;
1226: return 0;
1227: }
1229: /*@C
1230: DMMoabCreateElement - Adds an element of specified type to the primary set represented by the DM.
1232: Collective
1234: Input Parameters:
1235: + dm - The DM object
1236: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1237: . conn - The connectivity of the element
1238: - nverts - The number of vertices that form the element
1240: Output Parameter:
1241: . oelem - The handle to the element created and added to the DM object
1243: Level: beginner
1245: .seealso: DMMoabCreateSubmesh(), DMMoabCreateVertices()
1246: @*/
1247: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle* conn, PetscInt nverts, moab::EntityHandle* oelem)
1248: {
1249: moab::ErrorCode merr;
1250: DM_Moab *dmmoab;
1251: moab::EntityHandle elem;
1256: dmmoab = (DM_Moab*) dm->data;
1258: /* Insert new element */
1259: merr = dmmoab->mbiface->create_element(type, conn, nverts, elem); MBERRNM(merr);
1260: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1); MBERRNM(merr);
1262: if (oelem) *oelem = elem;
1263: return 0;
1264: }
1266: /*@C
1267: DMMoabCreateSubmesh - Creates a sub-DM object with a set that contains all vertices/elements of the parent
1268: in addition to providing support for dynamic mesh modifications. This is useful for AMR calculations to
1269: create a DM object on a refined level.
1271: Collective
1273: Input Parameters:
1274: . dm - The DM object
1276: Output Parameter:
1277: . newdm - The sub DM object with updated set information
1279: Level: advanced
1281: .seealso: DMCreate(), DMMoabCreateVertices(), DMMoabCreateElement()
1282: @*/
1283: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1284: {
1285: DM_Moab *dmmoab;
1286: DM_Moab *ndmmoab;
1287: moab::ErrorCode merr;
1291: dmmoab = (DM_Moab*) dm->data;
1293: /* Create the basic DMMoab object and keep the default parameters created by DM impls */
1294: DMMoabCreateMoab(((PetscObject)dm)->comm, dmmoab->mbiface, &dmmoab->ltog_tag, PETSC_NULL, newdm);
1296: /* get all the necessary handles from the private DM object */
1297: ndmmoab = (DM_Moab*) (*newdm)->data;
1299: /* set the sub-mesh's parent DM reference */
1300: ndmmoab->parent = &dm;
1302: /* create a file set to associate all entities in current mesh */
1303: merr = ndmmoab->mbiface->create_meshset(moab::MESHSET_SET, ndmmoab->fileset); MBERR("Creating file set failed", merr);
1305: /* create a meshset and then add old fileset as child */
1306: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal); MBERR("Adding child vertices to parent failed", merr);
1307: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal); MBERR("Adding child elements to parent failed", merr);
1309: /* preserve the field association between the parent and sub-mesh objects */
1310: DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames);
1311: return 0;
1312: }
1314: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1315: {
1316: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
1317: const char *name;
1318: MPI_Comm comm;
1319: PetscMPIInt size;
1321: PetscObjectGetComm((PetscObject)dm, &comm);
1322: MPI_Comm_size(comm, &size);
1323: PetscObjectGetName((PetscObject) dm, &name);
1324: PetscViewerASCIIPushTab(viewer);
1325: if (name) PetscViewerASCIIPrintf(viewer, "%s in %D dimensions:\n", name, dmmoab->dim);
1326: else PetscViewerASCIIPrintf(viewer, "Mesh in %D dimensions:\n", dmmoab->dim);
1327: /* print details about the global mesh */
1328: {
1329: PetscViewerASCIIPushTab(viewer);
1330: PetscViewerASCIIPrintf(viewer, "Sizes: cells=%D, vertices=%D, blocks=%D\n", dmmoab->nele, dmmoab->n, dmmoab->bs);
1331: /* print boundary data */
1332: PetscViewerASCIIPrintf(viewer, "Boundary trace:\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size());
1333: {
1334: PetscViewerASCIIPushTab(viewer);
1335: PetscViewerASCIIPrintf(viewer, "cells=%D, faces=%D, vertices=%D\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size());
1336: PetscViewerASCIIPopTab(viewer);
1337: }
1338: /* print field data */
1339: PetscViewerASCIIPrintf(viewer, "Fields: %D components\n", dmmoab->numFields);
1340: {
1341: PetscViewerASCIIPushTab(viewer);
1342: for (int i = 0; i < dmmoab->numFields; ++i) {
1343: PetscViewerASCIIPrintf(viewer, "[%D] - %s\n", i, dmmoab->fieldNames[i]);
1344: }
1345: PetscViewerASCIIPopTab(viewer);
1346: }
1347: PetscViewerASCIIPopTab(viewer);
1348: }
1349: PetscViewerASCIIPopTab(viewer);
1350: PetscViewerFlush(viewer);
1351: return 0;
1352: }
1354: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1355: {
1356: return 0;
1357: }
1359: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1360: {
1361: return 0;
1362: }
1364: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1365: {
1366: PetscBool iascii, ishdf5, isvtk;
1370: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
1371: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERVTK, &isvtk);
1372: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
1373: if (iascii) {
1374: DMMoabView_Ascii(dm, viewer);
1375: } else if (ishdf5) {
1376: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1377: PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ);
1378: DMMoabView_HDF5(dm, viewer);
1379: PetscViewerPopFormat(viewer);
1380: #else
1381: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1382: #endif
1383: }
1384: else if (isvtk) {
1385: DMMoabView_VTK(dm, viewer);
1386: }
1387: return 0;
1388: }
1390: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
1391: {
1392: dm->ops->view = DMView_Moab;
1393: dm->ops->load = NULL /* DMLoad_Moab */;
1394: dm->ops->setfromoptions = DMSetFromOptions_Moab;
1395: dm->ops->clone = DMClone_Moab;
1396: dm->ops->setup = DMSetUp_Moab;
1397: dm->ops->createlocalsection = NULL;
1398: dm->ops->createdefaultconstraints = NULL;
1399: dm->ops->createglobalvector = DMCreateGlobalVector_Moab;
1400: dm->ops->createlocalvector = DMCreateLocalVector_Moab;
1401: dm->ops->getlocaltoglobalmapping = NULL;
1402: dm->ops->createfieldis = NULL;
1403: dm->ops->createcoordinatedm = NULL /* DMCreateCoordinateDM_Moab */;
1404: dm->ops->getcoloring = NULL;
1405: dm->ops->creatematrix = DMCreateMatrix_Moab;
1406: dm->ops->createinterpolation = DMCreateInterpolation_Moab;
1407: dm->ops->createinjection = NULL /* DMCreateInjection_Moab */;
1408: dm->ops->refine = DMRefine_Moab;
1409: dm->ops->coarsen = DMCoarsen_Moab;
1410: dm->ops->refinehierarchy = DMRefineHierarchy_Moab;
1411: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Moab;
1412: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Moab;
1413: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Moab;
1414: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Moab;
1415: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Moab;
1416: dm->ops->destroy = DMDestroy_Moab;
1417: dm->ops->createsubdm = NULL /* DMCreateSubDM_Moab */;
1418: dm->ops->getdimpoints = NULL /* DMGetDimPoints_Moab */;
1419: dm->ops->locatepoints = NULL /* DMLocatePoints_Moab */;
1420: return 0;
1421: }
1423: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1424: {
1425: /* get all the necessary handles from the private DM object */
1426: (*newdm)->data = (DM_Moab*) dm->data;
1427: ((DM_Moab*)dm->data)->refct++;
1429: PetscObjectChangeTypeName((PetscObject) *newdm, DMMOAB);
1430: DMInitialize_Moab(*newdm);
1431: return 0;
1432: }
1434: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1435: {
1437: PetscNewLog(dm, (DM_Moab**)&dm->data);
1439: ((DM_Moab*)dm->data)->bs = 1;
1440: ((DM_Moab*)dm->data)->numFields = 1;
1441: ((DM_Moab*)dm->data)->n = 0;
1442: ((DM_Moab*)dm->data)->nloc = 0;
1443: ((DM_Moab*)dm->data)->nghost = 0;
1444: ((DM_Moab*)dm->data)->nele = 0;
1445: ((DM_Moab*)dm->data)->neleloc = 0;
1446: ((DM_Moab*)dm->data)->neleghost = 0;
1447: ((DM_Moab*)dm->data)->ltog_map = NULL;
1448: ((DM_Moab*)dm->data)->ltog_sendrecv = NULL;
1450: ((DM_Moab*)dm->data)->refct = 1;
1451: ((DM_Moab*)dm->data)->parent = NULL;
1452: ((DM_Moab*)dm->data)->vlocal = new moab::Range();
1453: ((DM_Moab*)dm->data)->vowned = new moab::Range();
1454: ((DM_Moab*)dm->data)->vghost = new moab::Range();
1455: ((DM_Moab*)dm->data)->elocal = new moab::Range();
1456: ((DM_Moab*)dm->data)->eghost = new moab::Range();
1458: DMInitialize_Moab(dm);
1459: return 0;
1460: }