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: Level: intermediate
19: Reference:
20: . * - https://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html
22: .seealso: `DMMOAB`, `DMType`, `DMMoabCreate()`, `DMCreate()`, `DMSetType()`, `DMMoabCreateMoab()`
23: M*/
25: /* External function declarations here */
26: PETSC_EXTERN PetscErrorCode DMCreateInterpolation_Moab(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
27: PETSC_EXTERN PetscErrorCode DMCreateDefaultConstraints_Moab(DM dm);
28: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
29: PETSC_EXTERN PetscErrorCode DMCreateCoordinateDM_Moab(DM dm, DM *cdm);
30: PETSC_EXTERN PetscErrorCode DMRefine_Moab(DM dm, MPI_Comm comm, DM *dmRefined);
31: PETSC_EXTERN PetscErrorCode DMCoarsen_Moab(DM dm, MPI_Comm comm, DM *dmCoarsened);
32: PETSC_EXTERN PetscErrorCode DMRefineHierarchy_Moab(DM dm, PetscInt nlevels, DM dmRefined[]);
33: PETSC_EXTERN PetscErrorCode DMCoarsenHierarchy_Moab(DM dm, PetscInt nlevels, DM dmCoarsened[]);
34: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm);
35: PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM, Vec *);
36: PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM, Vec *);
37: PETSC_EXTERN PetscErrorCode DMCreateMatrix_Moab(DM dm, Mat *J);
38: PETSC_EXTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM, Vec, InsertMode, Vec);
39: PETSC_EXTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM, Vec, InsertMode, Vec);
40: PETSC_EXTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM, Vec, InsertMode, Vec);
41: PETSC_EXTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM, Vec, InsertMode, Vec);
43: /* Un-implemented routines */
44: /*
45: PETSC_EXTERN PetscErrorCode DMCreatelocalsection_Moab(DM dm);
46: PETSC_EXTERN PetscErrorCode DMCreateInjection_Moab(DM dmCoarse, DM dmFine, Mat *mat);
47: PETSC_EXTERN PetscErrorCode DMLoad_Moab(DM dm, PetscViewer viewer);
48: PETSC_EXTERN PetscErrorCode DMGetDimPoints_Moab(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd);
49: PETSC_EXTERN PetscErrorCode DMCreateSubDM_Moab(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm);
50: PETSC_EXTERN PetscErrorCode DMLocatePoints_Moab(DM dm, Vec v, IS *cellIS);
51: */
53: /*@C
54: DMMoabCreate - Creates a `DMMOAB` object, which encapsulates a moab instance
56: Collective
58: Input Parameter:
59: . comm - The communicator for the `DMMOAB` object
61: Output Parameter:
62: . dmb - The `DMMOAB` object
64: Level: beginner
66: .seealso: `DMMOAB`, `DMMoabCreateMoab()`
67: @*/
68: PetscErrorCode DMMoabCreate(MPI_Comm comm, DM *dmb)
69: {
70: PetscFunctionBegin;
71: PetscAssertPointer(dmb, 2);
72: PetscCall(DMCreate(comm, dmb));
73: PetscCall(DMSetType(*dmb, DMMOAB));
74: PetscFunctionReturn(PETSC_SUCCESS);
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 `DMMOAB` 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: .seealso: `DMMOAB`, `DMMoabCreate()`
95: @*/
96: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
97: {
98: moab::ErrorCode merr;
99: DM dmmb;
100: DM_Moab *dmmoab;
102: PetscFunctionBegin;
103: PetscAssertPointer(dmb, 6);
105: PetscCall(DMMoabCreate(comm, &dmmb));
106: dmmoab = (DM_Moab *)(dmmb)->data;
108: if (!mbiface) {
109: dmmoab->mbiface = new moab::Core();
110: dmmoab->icreatedinstance = PETSC_TRUE;
111: } else {
112: dmmoab->mbiface = mbiface;
113: dmmoab->icreatedinstance = PETSC_FALSE;
114: }
116: /* by default the fileset = root set. This set stores the hierarchy of entities belonging to current DM */
117: dmmoab->fileset = 0;
118: dmmoab->hlevel = 0;
119: dmmoab->nghostrings = 0;
121: #ifdef MOAB_HAVE_MPI
122: moab::EntityHandle partnset;
124: /* Create root sets for each mesh. Then pass these
125: to the load_file functions to be populated. */
126: merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, partnset);
127: 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: PetscCall(PetscMalloc(dmmoab->numFields * sizeof(char *), &dmmoab->fieldNames));
137: PetscCall(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: PetscCall(DMMoabSetLocalToGlobalTag(dmmb, *ltog_tag));
148: } else {
149: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
150: 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);
155: MBERRNM(merr);
157: /* set the local range of entities (vertices) of interest */
158: if (range) PetscCall(DMMoabSetLocalVertices(dmmb, range));
159: *dmb = dmmb;
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: #ifdef MOAB_HAVE_MPI
165: /*@C
166: DMMoabGetParallelComm - Get the ParallelComm used with this `DMMOAB`
168: Collective
170: Input Parameter:
171: . dm - The `DMMOAB` object being set
173: Output Parameter:
174: . pcomm - The ParallelComm for the `DMMOAB`
176: Level: beginner
178: .seealso: `DMMOAB`, `DMMoabSetInterface()`
179: @*/
180: PetscErrorCode DMMoabGetParallelComm(DM dm, moab::ParallelComm **pcomm)
181: {
182: PetscFunctionBegin;
184: *pcomm = ((DM_Moab *)(dm)->data)->pcomm;
185: PetscFunctionReturn(PETSC_SUCCESS);
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: .seealso: `DMMOAB`, `DMMoabGetInterface()`
202: @*/
203: PetscErrorCode DMMoabSetInterface(DM dm, moab::Interface *mbiface)
204: {
205: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
207: PetscFunctionBegin;
209: PetscAssertPointer(mbiface, 2);
210: #ifdef MOAB_HAVE_MPI
211: dmmoab->pcomm = NULL;
212: #endif
213: dmmoab->mbiface = mbiface;
214: dmmoab->icreatedinstance = PETSC_FALSE;
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: /*@C
219: DMMoabGetInterface - Get the MOAB instance used with this `DMMOAB`
221: Collective
223: Input Parameter:
224: . dm - The `DMMOAB` object being set
226: Output Parameter:
227: . mbiface - The MOAB instance set on this `DMMOAB`
229: Level: beginner
231: .seealso: `DMMOAB`, `DMMoabSetInterface()`
232: @*/
233: PetscErrorCode DMMoabGetInterface(DM dm, moab::Interface **mbiface)
234: {
235: static PetscBool cite = PETSC_FALSE;
237: PetscFunctionBegin;
239: PetscCall(PetscCitationsRegister("@techreport{tautges_moab:_2004,\n type = {{SAND2004-1592}},\n title = {{MOAB:} A Mesh-Oriented Database}, institution = {Sandia National Laboratories},\n author = {Tautges, T. J. and Meyers, R. and Merkley, "
240: "K. and Stimpson, C. and Ernst, C.},\n year = {2004}, note = {Report}\n}\n",
241: &cite));
242: *mbiface = ((DM_Moab *)dm->data)->mbiface;
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
246: /*@C
247: DMMoabSetLocalVertices - Set the entities having DOFs on this `DMMOAB`
249: Collective
251: Input Parameters:
252: + dm - The `DMMOAB` object being set
253: - range - The entities treated by this `DMMOAB`
255: Level: beginner
257: .seealso: `DMMOAB`, `DMMoabGetAllVertices()`
258: @*/
259: PetscErrorCode DMMoabSetLocalVertices(DM dm, moab::Range *range)
260: {
261: moab::Range tmpvtxs;
262: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
264: PetscFunctionBegin;
266: dmmoab->vlocal->clear();
267: dmmoab->vowned->clear();
269: dmmoab->vlocal->insert(range->begin(), range->end());
271: #ifdef MOAB_HAVE_MPI
272: moab::ErrorCode merr;
273: /* filter based on parallel status */
274: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
275: MBERRNM(merr);
277: /* filter all the non-owned and shared entities out of the list */
278: tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
279: merr = dmmoab->pcomm->filter_pstatus(tmpvtxs, PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
280: MBERRNM(merr);
281: tmpvtxs = moab::subtract(tmpvtxs, *dmmoab->vghost);
282: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
283: #else
284: *dmmoab->vowned = *dmmoab->vlocal;
285: #endif
287: /* compute and cache the sizes of local and ghosted entities */
288: dmmoab->nloc = dmmoab->vowned->size();
289: dmmoab->nghost = dmmoab->vghost->size();
290: #ifdef MOAB_HAVE_MPI
291: PetscCall(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
292: #else
293: dmmoab->n = dmmoab->nloc;
294: #endif
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: /*@C
299: DMMoabGetAllVertices - Get the entities having DOFs on this `DMMOAB`
301: Collective
303: Input Parameter:
304: . dm - The `DMMOAB` object being set
306: Output Parameter:
307: . local - The local vertex entities in this `DMMOAB` = (owned+ghosted)
309: Level: beginner
311: .seealso: `DMMOAB`, `DMMoabGetLocalVertices()`
312: @*/
313: PetscErrorCode DMMoabGetAllVertices(DM dm, moab::Range *local)
314: {
315: PetscFunctionBegin;
317: if (local) *local = *((DM_Moab *)dm->data)->vlocal;
318: PetscFunctionReturn(PETSC_SUCCESS);
319: }
321: /*@C
322: DMMoabGetLocalVertices - Get the entities having DOFs on this `DMMOAB`
324: Collective
326: Input Parameter:
327: . dm - The `DMMOAB` object being set
329: Output Parameters:
330: + owned - The owned vertex entities in this `DMMOAB`
331: - ghost - The ghosted entities (non-owned) stored locally in this partition
333: Level: beginner
335: .seealso: `DMMOAB`, `DMMoabGetAllVertices()`
336: @*/
337: PetscErrorCode DMMoabGetLocalVertices(DM dm, const moab::Range **owned, const moab::Range **ghost)
338: {
339: PetscFunctionBegin;
341: if (owned) *owned = ((DM_Moab *)dm->data)->vowned;
342: if (ghost) *ghost = ((DM_Moab *)dm->data)->vghost;
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: /*@C
347: DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned
349: Collective
351: Input Parameter:
352: . dm - The `DMMOAB` object being set
354: Output Parameter:
355: . range - The entities owned locally
357: Level: beginner
359: .seealso: `DMMOAB`, `DMMoabSetLocalElements()`
360: @*/
361: PetscErrorCode DMMoabGetLocalElements(DM dm, const moab::Range **range)
362: {
363: PetscFunctionBegin;
365: if (range) *range = ((DM_Moab *)dm->data)->elocal;
366: PetscFunctionReturn(PETSC_SUCCESS);
367: }
369: /*@C
370: DMMoabSetLocalElements - Set the entities having DOFs on this `DMMOAB`
372: Collective
374: Input Parameters:
375: + dm - The `DMMOAB` object being set
376: - range - The entities treated by this `DMMOAB`
378: Level: beginner
380: .seealso: `DMMOAB`, `DMMoabGetLocalElements()`
381: @*/
382: PetscErrorCode DMMoabSetLocalElements(DM dm, moab::Range *range)
383: {
384: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
386: PetscFunctionBegin;
388: dmmoab->elocal->clear();
389: dmmoab->eghost->clear();
390: dmmoab->elocal->insert(range->begin(), range->end());
391: #ifdef MOAB_HAVE_MPI
392: moab::ErrorCode merr;
393: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
394: MBERRNM(merr);
395: *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
396: #endif
397: dmmoab->neleloc = dmmoab->elocal->size();
398: dmmoab->neleghost = dmmoab->eghost->size();
399: #ifdef MOAB_HAVE_MPI
400: PetscCall(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
401: PetscCall(PetscInfo(dm, "Created %" PetscInt_FMT " local and %" PetscInt_FMT " global elements.\n", dmmoab->neleloc, dmmoab->nele));
402: #else
403: dmmoab->nele = dmmoab->neleloc;
404: #endif
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: /*@C
409: DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
411: Collective
413: Input Parameters:
414: + dm - The `DMMOAB` object being set
415: - ltogtag - The `DMMOAB` tag used for local to global ids
417: Level: beginner
419: .seealso: `DMMOAB`, `DMMoabGetLocalToGlobalTag()`
420: @*/
421: PetscErrorCode DMMoabSetLocalToGlobalTag(DM dm, moab::Tag ltogtag)
422: {
423: PetscFunctionBegin;
425: ((DM_Moab *)dm->data)->ltog_tag = ltogtag;
426: PetscFunctionReturn(PETSC_SUCCESS);
427: }
429: /*@C
430: DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
432: Collective
434: Input Parameter:
435: . dm - The `DMMOAB` object being set
437: Output Parameter:
438: . ltog_tag - The MOAB tag used for local to global ids
440: Level: beginner
442: .seealso: `DMMOAB`, `DMMoabSetLocalToGlobalTag()`
443: @*/
444: PetscErrorCode DMMoabGetLocalToGlobalTag(DM dm, moab::Tag *ltog_tag)
445: {
446: PetscFunctionBegin;
448: *ltog_tag = ((DM_Moab *)dm->data)->ltog_tag;
449: PetscFunctionReturn(PETSC_SUCCESS);
450: }
452: /*@C
453: DMMoabSetBlockSize - Set the block size used with this `DMMOAB`
455: Collective
457: Input Parameters:
458: + dm - The `DMMOAB` object being set
459: - bs - The block size used with this `DMMOAB`
461: Level: beginner
463: .seealso: `DMMOAB`, `DMMoabGetBlockSize()`
464: @*/
465: PetscErrorCode DMMoabSetBlockSize(DM dm, PetscInt bs)
466: {
467: PetscFunctionBegin;
469: ((DM_Moab *)dm->data)->bs = bs;
470: PetscFunctionReturn(PETSC_SUCCESS);
471: }
473: /*@C
474: DMMoabGetBlockSize - Get the block size used with this `DMMOAB`
476: Collective
478: Input Parameter:
479: . dm - The `DMMOAB` object being set
481: Output Parameter:
482: . bs - The block size used with this `DMMOAB`
484: Level: beginner
486: .seealso: `DMMOAB`, `DMMoabSetBlockSize()`
487: @*/
488: PetscErrorCode DMMoabGetBlockSize(DM dm, PetscInt *bs)
489: {
490: PetscFunctionBegin;
492: *bs = ((DM_Moab *)dm->data)->bs;
493: PetscFunctionReturn(PETSC_SUCCESS);
494: }
496: /*@C
497: DMMoabGetSize - Get the global vertex size used with this `DMMOAB`
499: Collective
501: Input Parameter:
502: . dm - The `DMMOAB` object being set
504: Output Parameters:
505: + neg - The number of global elements in the `DMMOAB` instance
506: - nvg - The number of global vertices in the `DMMOAB` instance
508: Level: beginner
510: .seealso: `DMMOAB`, `DMMoabGetLocalSize()`
511: @*/
512: PetscErrorCode DMMoabGetSize(DM dm, PetscInt *neg, PetscInt *nvg)
513: {
514: PetscFunctionBegin;
516: if (neg) *neg = ((DM_Moab *)dm->data)->nele;
517: if (nvg) *nvg = ((DM_Moab *)dm->data)->n;
518: PetscFunctionReturn(PETSC_SUCCESS);
519: }
521: /*@C
522: DMMoabGetLocalSize - Get the local and ghosted vertex size used with this `DMMOAB`
524: Collective
526: Input Parameter:
527: . dm - The `DMMOAB` object being set
529: Output Parameters:
530: + nel - The number of owned elements in this processor
531: . neg - The number of ghosted elements in this processor
532: . nvl - The number of owned vertices in this processor
533: - nvg - The number of ghosted vertices in this processor
535: Level: beginner
537: .seealso: `DMMOAB`, `DMMoabGetSize()`
538: @*/
539: PetscErrorCode DMMoabGetLocalSize(DM dm, PetscInt *nel, PetscInt *neg, PetscInt *nvl, PetscInt *nvg)
540: {
541: PetscFunctionBegin;
543: if (nel) *nel = ((DM_Moab *)dm->data)->neleloc;
544: if (neg) *neg = ((DM_Moab *)dm->data)->neleghost;
545: if (nvl) *nvl = ((DM_Moab *)dm->data)->nloc;
546: if (nvg) *nvg = ((DM_Moab *)dm->data)->nghost;
547: PetscFunctionReturn(PETSC_SUCCESS);
548: }
550: /*@C
551: DMMoabGetOffset - Get the local offset for the global vector
553: Collective
555: Input Parameter:
556: . dm - The `DMMOAB` object being set
558: Output Parameter:
559: . offset - The local offset for the global vector
561: Level: beginner
563: .seealso: `DMMOAB`, `DMMoabGetDimension()`
564: @*/
565: PetscErrorCode DMMoabGetOffset(DM dm, PetscInt *offset)
566: {
567: PetscFunctionBegin;
569: *offset = ((DM_Moab *)dm->data)->vstart;
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: /*@C
574: DMMoabGetDimension - Get the dimension of the `DM` Mesh
576: Collective
578: Input Parameter:
579: . dm - The `DMMOAB` object
581: Output Parameter:
582: . dim - The dimension of `DM`
584: Level: beginner
586: .seealso: `DMMOAB`, `DMMoabGetOffset()`
587: @*/
588: PetscErrorCode DMMoabGetDimension(DM dm, PetscInt *dim)
589: {
590: PetscFunctionBegin;
592: *dim = ((DM_Moab *)dm->data)->dim;
593: PetscFunctionReturn(PETSC_SUCCESS);
594: }
596: /*@C
597: DMMoabGetHierarchyLevel - Get the current level of the mesh hierarchy
598: generated through uniform refinement.
600: Collective
602: Input Parameter:
603: . dm - The `DMMOAB` object being set
605: Output Parameter:
606: . nlevel - The current mesh hierarchy level
608: Level: beginner
610: .seealso: `DMMOAB`, `DMMoabGetMaterialBlock()`
611: @*/
612: PetscErrorCode DMMoabGetHierarchyLevel(DM dm, PetscInt *nlevel)
613: {
614: PetscFunctionBegin;
616: if (nlevel) *nlevel = ((DM_Moab *)dm->data)->hlevel;
617: PetscFunctionReturn(PETSC_SUCCESS);
618: }
620: /*@C
621: DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh
623: Collective
625: Input Parameters:
626: + dm - The `DMMOAB` object
627: - ehandle - The element entity handle
629: Output Parameter:
630: . mat - The material ID for the current entity
632: Level: beginner
634: .seealso: `DMMOAB`, `DMMoabGetHierarchyLevel()`
635: @*/
636: PetscErrorCode DMMoabGetMaterialBlock(DM dm, const moab::EntityHandle ehandle, PetscInt *mat)
637: {
638: DM_Moab *dmmoab;
640: PetscFunctionBegin;
642: if (*mat) {
643: dmmoab = (DM_Moab *)(dm)->data;
644: *mat = dmmoab->materials[dmmoab->elocal->index(ehandle)];
645: }
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
649: /*@C
650: DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities
652: Collective
654: Input Parameters:
655: + dm - The `DMMOAB` object
656: . nconn - Number of entities whose coordinates are needed
657: - conn - The vertex entity handles
659: Output Parameter:
660: . vpos - The coordinates of the requested vertex entities
662: Level: beginner
664: .seealso: `DMMOAB`, `DMMoabGetVertexConnectivity()`
665: @*/
666: PetscErrorCode DMMoabGetVertexCoordinates(DM dm, PetscInt nconn, const moab::EntityHandle *conn, PetscReal *vpos)
667: {
668: DM_Moab *dmmoab;
669: moab::ErrorCode merr;
671: PetscFunctionBegin;
673: PetscAssertPointer(conn, 3);
674: PetscAssertPointer(vpos, 4);
675: dmmoab = (DM_Moab *)(dm)->data;
677: /* Get connectivity information in MOAB canonical ordering */
678: if (dmmoab->hlevel) {
679: merr = dmmoab->hierarchy->get_coordinates(const_cast<moab::EntityHandle *>(conn), nconn, dmmoab->hlevel, vpos);
680: MBERRNM(merr);
681: } else {
682: merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);
683: MBERRNM(merr);
684: }
685: PetscFunctionReturn(PETSC_SUCCESS);
686: }
688: /*@C
689: DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity
691: Collective
693: Input Parameters:
694: + dm - The `DMMOAB` object
695: - vhandle - Vertex entity handle
697: Output Parameters:
698: + nconn - Number of entities whose coordinates are needed
699: - conn - The vertex entity handles
701: Level: beginner
703: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabRestoreVertexConnectivity()`
704: @*/
705: PetscErrorCode DMMoabGetVertexConnectivity(DM dm, moab::EntityHandle vhandle, PetscInt *nconn, moab::EntityHandle **conn)
706: {
707: DM_Moab *dmmoab;
708: std::vector<moab::EntityHandle> adj_entities, connect;
709: moab::ErrorCode merr;
711: PetscFunctionBegin;
713: PetscAssertPointer(conn, 4);
714: dmmoab = (DM_Moab *)(dm)->data;
716: /* Get connectivity information in MOAB canonical ordering */
717: merr = dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION);
718: MBERRNM(merr);
719: merr = dmmoab->mbiface->get_connectivity(&adj_entities[0], adj_entities.size(), connect);
720: MBERRNM(merr);
722: if (conn) {
723: PetscCall(PetscMalloc(sizeof(moab::EntityHandle) * connect.size(), conn));
724: PetscCall(PetscArraycpy(*conn, &connect[0], connect.size()));
725: }
726: if (nconn) *nconn = connect.size();
727: PetscFunctionReturn(PETSC_SUCCESS);
728: }
730: /*@C
731: DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity
733: Collective
735: Input Parameters:
736: + dm - The `DMMOAB` object
737: . ehandle - Vertex entity handle
738: . nconn - Number of entities whose coordinates are needed
739: - conn - The vertex entity handles
741: Level: beginner
743: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`
744: @*/
745: PetscErrorCode DMMoabRestoreVertexConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, moab::EntityHandle **conn)
746: {
747: PetscFunctionBegin;
749: PetscAssertPointer(conn, 4);
751: if (conn) PetscCall(PetscFree(*conn));
752: if (nconn) *nconn = 0;
753: PetscFunctionReturn(PETSC_SUCCESS);
754: }
756: /*@C
757: DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity
759: Collective
761: Input Parameters:
762: + dm - The `DMMOAB` object
763: - ehandle - Vertex entity handle
765: Output Parameters:
766: + nconn - Number of entities whose coordinates are needed
767: - conn - The vertex entity handles
769: Level: beginner
771: .seealso: `DMMOAB`, `DMMoabGetVertexCoordinates()`, `DMMoabGetVertexConnectivity()`, `DMMoabRestoreVertexConnectivity()`
772: @*/
773: PetscErrorCode DMMoabGetElementConnectivity(DM dm, moab::EntityHandle ehandle, PetscInt *nconn, const moab::EntityHandle **conn)
774: {
775: DM_Moab *dmmoab;
776: const moab::EntityHandle *connect;
777: std::vector<moab::EntityHandle> vconn;
778: moab::ErrorCode merr;
779: PetscInt nnodes;
781: PetscFunctionBegin;
783: PetscAssertPointer(conn, 4);
784: dmmoab = (DM_Moab *)(dm)->data;
786: /* Get connectivity information in MOAB canonical ordering */
787: merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes);
788: MBERRNM(merr);
789: if (conn) *conn = connect;
790: if (nconn) *nconn = nnodes;
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }
794: /*@C
795: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
797: Collective
799: Input Parameters:
800: + dm - The `DMMOAB` object
801: - ent - Entity handle
803: Output Parameter:
804: . ent_on_boundary - `PETSC_TRUE` if entity on boundary; `PETSC_FALSE` otherwise
806: Level: beginner
808: .seealso: `DMMOAB`, `DMMoabCheckBoundaryVertices()`
809: @*/
810: PetscErrorCode DMMoabIsEntityOnBoundary(DM dm, const moab::EntityHandle ent, PetscBool *ent_on_boundary)
811: {
812: moab::EntityType etype;
813: DM_Moab *dmmoab;
814: PetscInt edim;
816: PetscFunctionBegin;
818: PetscAssertPointer(ent_on_boundary, 3);
819: dmmoab = (DM_Moab *)(dm)->data;
821: /* get the entity type and handle accordingly */
822: etype = dmmoab->mbiface->type_from_handle(ent);
823: PetscCheck(etype < moab::MBPOLYHEDRON, PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Entity type on the boundary skin is invalid. EntityType = %" PetscInt_FMT, etype);
825: /* get the entity dimension */
826: edim = dmmoab->mbiface->dimension_from_handle(ent);
828: *ent_on_boundary = PETSC_FALSE;
829: if (etype == moab::MBVERTEX && edim == 0) {
830: *ent_on_boundary = ((dmmoab->bndyvtx->index(ent) >= 0) ? PETSC_TRUE : PETSC_FALSE);
831: } else {
832: if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
833: if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
834: } else { /* next check the lower-dimensional faces */
835: if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary = PETSC_TRUE;
836: }
837: }
838: PetscFunctionReturn(PETSC_SUCCESS);
839: }
841: /*@C
842: DMMoabCheckBoundaryVertices - Check whether a given entity is on the boundary (vertex, edge, face, element)
844: Input Parameters:
845: + dm - The `DMMOAB` object
846: . nconn - Number of handles
847: - cnt - Array of entity handles
849: Output Parameter:
850: . isbdvtx - Array of boundary markers - `PETSC_TRUE` if entity on boundary; `PETSC_FALSE` otherwise
852: Level: beginner
854: .seealso: `DMMOAB`, `DMMoabIsEntityOnBoundary()`
855: @*/
856: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm, PetscInt nconn, const moab::EntityHandle *cnt, PetscBool *isbdvtx)
857: {
858: DM_Moab *dmmoab;
859: PetscInt i;
861: PetscFunctionBegin;
863: PetscAssertPointer(cnt, 3);
864: PetscAssertPointer(isbdvtx, 4);
865: dmmoab = (DM_Moab *)(dm)->data;
867: for (i = 0; i < nconn; ++i) isbdvtx[i] = (dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE : PETSC_FALSE);
868: PetscFunctionReturn(PETSC_SUCCESS);
869: }
871: /*@C
872: DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary
874: Input Parameter:
875: . dm - The `DMMOAB` object
877: Output Parameters:
878: + bdvtx - Boundary vertices
879: . bdelems - Boundary elements
880: - bdfaces - Boundary faces
882: Level: beginner
884: .seealso: `DMMOAB`, `DMMoabCheckBoundaryVertices()`, `DMMoabIsEntityOnBoundary()`
885: @*/
886: PetscErrorCode DMMoabGetBoundaryMarkers(DM dm, const moab::Range **bdvtx, const moab::Range **bdelems, const moab::Range **bdfaces)
887: {
888: DM_Moab *dmmoab;
890: PetscFunctionBegin;
892: dmmoab = (DM_Moab *)(dm)->data;
894: if (bdvtx) *bdvtx = dmmoab->bndyvtx;
895: if (bdfaces) *bdfaces = dmmoab->bndyfaces;
896: if (bdelems) *bdfaces = dmmoab->bndyelems;
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)
901: {
902: PetscInt i;
903: moab::ErrorCode merr;
904: DM_Moab *dmmoab = (DM_Moab *)dm->data;
906: PetscFunctionBegin;
909: dmmoab->refct--;
910: if (!dmmoab->refct) {
911: delete dmmoab->vlocal;
912: delete dmmoab->vowned;
913: delete dmmoab->vghost;
914: delete dmmoab->elocal;
915: delete dmmoab->eghost;
916: delete dmmoab->bndyvtx;
917: delete dmmoab->bndyfaces;
918: delete dmmoab->bndyelems;
920: PetscCall(PetscFree(dmmoab->gsindices));
921: PetscCall(PetscFree2(dmmoab->gidmap, dmmoab->lidmap));
922: PetscCall(PetscFree(dmmoab->dfill));
923: PetscCall(PetscFree(dmmoab->ofill));
924: PetscCall(PetscFree(dmmoab->materials));
925: if (dmmoab->fieldNames) {
926: for (i = 0; i < dmmoab->numFields; i++) PetscCall(PetscFree(dmmoab->fieldNames[i]));
927: PetscCall(PetscFree(dmmoab->fieldNames));
928: }
930: if (dmmoab->nhlevels) {
931: PetscCall(PetscFree(dmmoab->hsets));
932: dmmoab->nhlevels = 0;
933: if (!dmmoab->hlevel && dmmoab->icreatedinstance) delete dmmoab->hierarchy;
934: dmmoab->hierarchy = NULL;
935: }
937: if (dmmoab->icreatedinstance) {
938: delete dmmoab->pcomm;
939: merr = dmmoab->mbiface->delete_mesh();
940: MBERRNM(merr);
941: delete dmmoab->mbiface;
942: }
943: dmmoab->mbiface = NULL;
944: #ifdef MOAB_HAVE_MPI
945: dmmoab->pcomm = NULL;
946: #endif
947: PetscCall(VecScatterDestroy(&dmmoab->ltog_sendrecv));
948: PetscCall(ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map));
949: PetscCall(PetscFree(dm->data));
950: }
951: PetscFunctionReturn(PETSC_SUCCESS);
952: }
954: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(DM dm, PetscOptionItems *PetscOptionsObject)
955: {
956: DM_Moab *dmmoab = (DM_Moab *)dm->data;
958: PetscFunctionBegin;
959: PetscOptionsHeadBegin(PetscOptionsObject, "DMMoab Options");
960: PetscCall(PetscOptionsBoundedInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL, 0));
961: PetscCall(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));
962: /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
963: PetscCall(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));
964: PetscCall(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));
965: PetscCall(PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum *)&dmmoab->read_mode, NULL));
966: PetscCall(PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum *)&dmmoab->write_mode, NULL));
967: PetscOptionsHeadEnd();
968: PetscFunctionReturn(PETSC_SUCCESS);
969: }
971: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)
972: {
973: moab::ErrorCode merr;
974: Vec local, global;
975: IS from, to;
976: moab::Range::iterator iter;
977: PetscInt i, j, f, bs, vent, totsize, *lgmap;
978: DM_Moab *dmmoab = (DM_Moab *)dm->data;
979: moab::Range adjs;
981: PetscFunctionBegin;
983: /* Get the local and shared vertices and cache it */
984: PetscCheck(dmmoab->mbiface != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB Interface before calling SetUp.");
985: #ifdef MOAB_HAVE_MPI
986: PetscCheck(dmmoab->pcomm != NULL, PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB ParallelComm object before calling SetUp.");
987: #endif
989: /* Get the entities recursively in the current part of the mesh, if user did not set the local vertices explicitly */
990: if (dmmoab->vlocal->empty()) {
991: //merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,*dmmoab->vlocal,true);MBERRNM(merr);
992: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, 0, *dmmoab->vlocal, false);
993: MBERRNM(merr);
995: #ifdef MOAB_HAVE_MPI
996: /* filter based on parallel status */
997: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, dmmoab->vowned);
998: MBERRNM(merr);
1000: /* filter all the non-owned and shared entities out of the list */
1001: // *dmmoab->vghost = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
1002: adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
1003: merr = dmmoab->pcomm->filter_pstatus(adjs, PSTATUS_GHOST | PSTATUS_INTERFACE, PSTATUS_OR, -1, dmmoab->vghost);
1004: MBERRNM(merr);
1005: adjs = moab::subtract(adjs, *dmmoab->vghost);
1006: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
1007: #else
1008: *dmmoab->vowned = *dmmoab->vlocal;
1009: #endif
1011: /* compute and cache the sizes of local and ghosted entities */
1012: dmmoab->nloc = dmmoab->vowned->size();
1013: dmmoab->nghost = dmmoab->vghost->size();
1015: #ifdef MOAB_HAVE_MPI
1016: PetscCall(MPIU_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
1017: PetscCall(PetscInfo(NULL, "Filset ID: %lu, Vertices: local - %zu, owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->fileset, dmmoab->vlocal->size(), dmmoab->nloc, dmmoab->nghost));
1018: #else
1019: dmmoab->n = dmmoab->nloc;
1020: #endif
1021: }
1023: {
1024: /* get the information about the local elements in the mesh */
1025: dmmoab->eghost->clear();
1027: /* first decipher the leading dimension */
1028: for (i = 3; i > 0; i--) {
1029: dmmoab->elocal->clear();
1030: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, false);
1031: MBERRNM(merr);
1033: /* store the current mesh dimension */
1034: if (dmmoab->elocal->size()) {
1035: dmmoab->dim = i;
1036: break;
1037: }
1038: }
1040: PetscCall(DMSetDimension(dm, dmmoab->dim));
1042: #ifdef MOAB_HAVE_MPI
1043: /* filter the ghosted and owned element list */
1044: *dmmoab->eghost = *dmmoab->elocal;
1045: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1046: MBERRNM(merr);
1047: *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
1048: #endif
1050: dmmoab->neleloc = dmmoab->elocal->size();
1051: dmmoab->neleghost = dmmoab->eghost->size();
1053: #ifdef MOAB_HAVE_MPI
1054: PetscCall(MPIU_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm));
1055: PetscCall(PetscInfo(NULL, "%d-dim elements: owned - %" PetscInt_FMT ", ghosted - %" PetscInt_FMT ".\n", dmmoab->dim, dmmoab->neleloc, dmmoab->neleghost));
1056: #else
1057: dmmoab->nele = dmmoab->neleloc;
1058: #endif
1059: }
1061: bs = dmmoab->bs;
1062: if (!dmmoab->ltog_tag) {
1063: /* Get the global ID tag. The global ID tag is applied to each
1064: vertex. It acts as an global identifier which MOAB uses to
1065: assemble the individual pieces of the mesh */
1066: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);
1067: MBERRNM(merr);
1068: }
1070: totsize = dmmoab->vlocal->size();
1071: PetscCheck(totsize == dmmoab->nloc + dmmoab->nghost, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Mismatch between local and owned+ghost vertices. %" PetscInt_FMT " != %" PetscInt_FMT ".", totsize, dmmoab->nloc + dmmoab->nghost);
1072: PetscCall(PetscCalloc1(totsize, &dmmoab->gsindices));
1073: {
1074: /* first get the local indices */
1075: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vowned, &dmmoab->gsindices[0]);
1076: MBERRNM(merr);
1077: if (dmmoab->nghost) { /* next get the ghosted indices */
1078: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag, *dmmoab->vghost, &dmmoab->gsindices[dmmoab->nloc]);
1079: MBERRNM(merr);
1080: }
1082: /* find out the local and global minima of GLOBAL_ID */
1083: dmmoab->lminmax[0] = dmmoab->lminmax[1] = dmmoab->gsindices[0];
1084: for (i = 0; i < totsize; ++i) {
1085: if (dmmoab->lminmax[0] > dmmoab->gsindices[i]) dmmoab->lminmax[0] = dmmoab->gsindices[i];
1086: if (dmmoab->lminmax[1] < dmmoab->gsindices[i]) dmmoab->lminmax[1] = dmmoab->gsindices[i];
1087: }
1089: PetscCall(MPIU_Allreduce(&dmmoab->lminmax[0], &dmmoab->gminmax[0], 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm));
1090: PetscCall(MPIU_Allreduce(&dmmoab->lminmax[1], &dmmoab->gminmax[1], 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm));
1092: /* set the GID map */
1093: for (i = 0; i < totsize; ++i) { dmmoab->gsindices[i] -= dmmoab->gminmax[0]; /* zero based index needed for IS */ }
1094: dmmoab->lminmax[0] -= dmmoab->gminmax[0];
1095: dmmoab->lminmax[1] -= dmmoab->gminmax[0];
1097: PetscCall(PetscInfo(NULL, "GLOBAL_ID: Local [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "], Global [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "]\n", dmmoab->lminmax[0], dmmoab->lminmax[1], dmmoab->gminmax[0], dmmoab->gminmax[1]));
1098: }
1099: PetscCheck(dmmoab->bs == dmmoab->numFields || dmmoab->bs == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mismatch between block size and number of component fields. %" PetscInt_FMT " != 1 OR %" PetscInt_FMT " != %" PetscInt_FMT ".", dmmoab->bs, dmmoab->bs,
1100: dmmoab->numFields);
1102: {
1103: dmmoab->seqstart = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->front());
1104: dmmoab->seqend = dmmoab->mbiface->id_from_handle(dmmoab->vlocal->back());
1105: PetscCall(PetscInfo(NULL, "SEQUENCE: Local [min, max] - [%" PetscInt_FMT ", %" PetscInt_FMT "]\n", dmmoab->seqstart, dmmoab->seqend));
1107: PetscCall(PetscMalloc2(dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->gidmap, dmmoab->seqend - dmmoab->seqstart + 1, &dmmoab->lidmap));
1108: PetscCall(PetscMalloc1(totsize * dmmoab->numFields, &lgmap));
1110: i = j = 0;
1111: /* set the owned vertex data first */
1112: for (moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++, i++) {
1113: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1114: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1115: dmmoab->lidmap[vent] = i;
1116: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1117: }
1118: /* next arrange all the ghosted data information */
1119: for (moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++, i++) {
1120: vent = dmmoab->mbiface->id_from_handle(*iter) - dmmoab->seqstart;
1121: dmmoab->gidmap[vent] = dmmoab->gsindices[i];
1122: dmmoab->lidmap[vent] = i;
1123: for (f = 0; f < dmmoab->numFields; f++, j++) lgmap[j] = (bs > 1 ? dmmoab->gsindices[i] * dmmoab->numFields + f : totsize * f + dmmoab->gsindices[i]);
1124: }
1126: /* We need to create the Global to Local Vector Scatter Contexts
1127: 1) First create a local and global vector
1128: 2) Create a local and global IS
1129: 3) Create VecScatter and LtoGMapping objects
1130: 4) Cleanup the IS and Vec objects
1131: */
1132: PetscCall(DMCreateGlobalVector(dm, &global));
1133: PetscCall(DMCreateLocalVector(dm, &local));
1135: PetscCall(VecGetOwnershipRange(global, &dmmoab->vstart, &dmmoab->vend));
1137: /* global to local must retrieve ghost points */
1138: PetscCall(ISCreateStride(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, dmmoab->vstart, 1, &from));
1139: PetscCall(ISSetBlockSize(from, bs));
1141: PetscCall(ISCreateGeneral(((PetscObject)dm)->comm, dmmoab->nloc * dmmoab->numFields, &lgmap[0], PETSC_COPY_VALUES, &to));
1142: PetscCall(ISSetBlockSize(to, bs));
1144: if (!dmmoab->ltog_map) {
1145: /* create to the local to global mapping for vectors in order to use VecSetValuesLocal */
1146: PetscCall(ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm, dmmoab->bs, totsize * dmmoab->numFields, lgmap, PETSC_COPY_VALUES, &dmmoab->ltog_map));
1147: }
1149: /* now create the scatter object from local to global vector */
1150: PetscCall(VecScatterCreate(local, from, global, to, &dmmoab->ltog_sendrecv));
1152: /* clean up IS, Vec */
1153: PetscCall(PetscFree(lgmap));
1154: PetscCall(ISDestroy(&from));
1155: PetscCall(ISDestroy(&to));
1156: PetscCall(VecDestroy(&local));
1157: PetscCall(VecDestroy(&global));
1158: }
1160: dmmoab->bndyvtx = new moab::Range();
1161: dmmoab->bndyfaces = new moab::Range();
1162: dmmoab->bndyelems = new moab::Range();
1163: /* skin the boundary and store nodes */
1164: if (!dmmoab->hlevel) {
1165: /* get the skin vertices of boundary faces for the current partition and then filter
1166: the local, boundary faces, vertices and elements alone via PSTATUS flags;
1167: this should not give us any ghosted boundary, but if user needs such a functionality
1168: it would be easy to add it based on the find_skin query below */
1169: moab::Skinner skinner(dmmoab->mbiface);
1171: /* get the entities on the skin - only the faces */
1172: merr = skinner.find_skin(dmmoab->fileset, *dmmoab->elocal, false, *dmmoab->bndyfaces, NULL, true, true, false);
1173: MBERRNM(merr); // 'false' param indicates we want faces back, not vertices
1175: #ifdef MOAB_HAVE_MPI
1176: /* filter all the non-owned and shared entities out of the list */
1177: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1178: MBERRNM(merr);
1179: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_INTERFACE, PSTATUS_NOT);
1180: MBERRNM(merr);
1181: #endif
1183: /* get all the nodes via connectivity and the parent elements via adjacency information */
1184: merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false);
1185: MBERRNM(merr);
1186: merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyvtx, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION);
1187: MBERRNM(merr);
1188: } else {
1189: /* Let us query the hierarchy manager and get the results directly for this level */
1190: for (moab::Range::iterator iter = dmmoab->elocal->begin(); iter != dmmoab->elocal->end(); iter++) {
1191: moab::EntityHandle elemHandle = *iter;
1192: if (dmmoab->hierarchy->is_entity_on_boundary(elemHandle)) {
1193: dmmoab->bndyelems->insert(elemHandle);
1194: /* For this boundary element, query the vertices and add them to the list */
1195: std::vector<moab::EntityHandle> connect;
1196: merr = dmmoab->hierarchy->get_connectivity(elemHandle, dmmoab->hlevel, connect);
1197: MBERRNM(merr);
1198: for (unsigned iv = 0; iv < connect.size(); ++iv)
1199: if (dmmoab->hierarchy->is_entity_on_boundary(connect[iv])) dmmoab->bndyvtx->insert(connect[iv]);
1200: /* Next, let us query the boundary faces and add them also to the list */
1201: std::vector<moab::EntityHandle> faces;
1202: merr = dmmoab->hierarchy->get_adjacencies(elemHandle, dmmoab->dim - 1, faces);
1203: MBERRNM(merr);
1204: for (unsigned ifa = 0; ifa < faces.size(); ++ifa)
1205: if (dmmoab->hierarchy->is_entity_on_boundary(faces[ifa])) dmmoab->bndyfaces->insert(faces[ifa]);
1206: }
1207: }
1208: #ifdef MOAB_HAVE_MPI
1209: /* filter all the non-owned and shared entities out of the list */
1210: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyvtx, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1211: MBERRNM(merr);
1212: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1213: MBERRNM(merr);
1214: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1215: MBERRNM(merr);
1216: #endif
1217: }
1218: PetscCall(PetscInfo(NULL, "Found %zu boundary vertices, %zu boundary faces and %zu boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyfaces->size(), dmmoab->bndyelems->size()));
1220: /* Get the material sets and populate the data for all locally owned elements */
1221: {
1222: PetscCall(PetscCalloc1(dmmoab->elocal->size(), &dmmoab->materials));
1223: /* Get the count of entities of particular type from dmmoab->elocal
1224: -- Then, for each non-zero type, loop through and query the fileset to get the material tag data */
1225: moab::Range msets;
1226: merr = dmmoab->mbiface->get_entities_by_type_and_tag(dmmoab->fileset, moab::MBENTITYSET, &dmmoab->material_tag, NULL, 1, msets, moab::Interface::UNION);
1227: MBERRNM(merr);
1228: if (msets.size() == 0) PetscCall(PetscInfo(NULL, "No material sets found in the fileset.\n"));
1230: for (unsigned i = 0; i < msets.size(); ++i) {
1231: moab::Range msetelems;
1232: merr = dmmoab->mbiface->get_entities_by_dimension(msets[i], dmmoab->dim, msetelems, true);
1233: MBERRNM(merr);
1234: #ifdef MOAB_HAVE_MPI
1235: /* filter all the non-owned and shared entities out of the list */
1236: merr = dmmoab->pcomm->filter_pstatus(msetelems, PSTATUS_NOT_OWNED, PSTATUS_NOT);
1237: MBERRNM(merr);
1238: #endif
1240: int partID;
1241: moab::EntityHandle mset = msets[i];
1242: merr = dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &mset, 1, &partID);
1243: MBERRNM(merr);
1245: for (unsigned j = 0; j < msetelems.size(); ++j) dmmoab->materials[dmmoab->elocal->index(msetelems[j])] = partID;
1246: }
1247: }
1249: PetscFunctionReturn(PETSC_SUCCESS);
1250: }
1252: /*@C
1253: DMMoabCreateVertices - Creates and adds several vertices to the primary set represented by the DM.
1255: Collective
1257: Input Parameters:
1258: + dm - The `DM` object
1259: . coords - The connectivity of the element
1260: - nverts - The number of vertices that form the element
1262: Output Parameter:
1263: . overts - The list of vertices that were created (can be `NULL`)
1265: Level: beginner
1267: .seealso: `DMMOAB`, `DMMoabCreateSubmesh()`, `DMMoabCreateElement()`
1268: @*/
1269: PetscErrorCode DMMoabCreateVertices(DM dm, const PetscReal *coords, PetscInt nverts, moab::Range *overts)
1270: {
1271: moab::ErrorCode merr;
1272: DM_Moab *dmmoab;
1273: moab::Range verts;
1275: PetscFunctionBegin;
1277: PetscAssertPointer(coords, 2);
1279: dmmoab = (DM_Moab *)dm->data;
1281: /* Insert new points */
1282: merr = dmmoab->mbiface->create_vertices(&coords[0], nverts, verts);
1283: MBERRNM(merr);
1284: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, verts);
1285: MBERRNM(merr);
1287: if (overts) *overts = verts;
1288: PetscFunctionReturn(PETSC_SUCCESS);
1289: }
1291: /*@C
1292: DMMoabCreateElement - Adds an element of specified type to the primary set represented by the DM.
1294: Collective
1296: Input Parameters:
1297: + dm - The DM object
1298: . type - The type of element to create and add (Edge/Tri/Quad/Tet/Hex/Prism/Pyramid/Polygon/Polyhedra)
1299: . conn - The connectivity of the element
1300: - nverts - The number of vertices that form the element
1302: Output Parameter:
1303: . oelem - The handle to the element created and added to the `DM` object
1305: Level: beginner
1307: .seealso: `DMMOAB`, `DMMoabCreateSubmesh()`, `DMMoabCreateVertices()`
1308: @*/
1309: PetscErrorCode DMMoabCreateElement(DM dm, const moab::EntityType type, const moab::EntityHandle *conn, PetscInt nverts, moab::EntityHandle *oelem)
1310: {
1311: moab::ErrorCode merr;
1312: DM_Moab *dmmoab;
1313: moab::EntityHandle elem;
1315: PetscFunctionBegin;
1317: PetscAssertPointer(conn, 3);
1319: dmmoab = (DM_Moab *)dm->data;
1321: /* Insert new element */
1322: merr = dmmoab->mbiface->create_element(type, conn, nverts, elem);
1323: MBERRNM(merr);
1324: merr = dmmoab->mbiface->add_entities(dmmoab->fileset, &elem, 1);
1325: MBERRNM(merr);
1327: if (oelem) *oelem = elem;
1328: PetscFunctionReturn(PETSC_SUCCESS);
1329: }
1331: /*@C
1332: DMMoabCreateSubmesh - Creates a sub-`DM` object with a set that contains all vertices/elements of the parent
1333: in addition to providing support for dynamic mesh modifications. This is useful for AMR calculations to
1334: create a DM object on a refined level.
1336: Collective
1338: Input Parameters:
1339: . dm - The `DM` object
1341: Output Parameter:
1342: . newdm - The sub `DM` object with updated set information
1344: Level: advanced
1346: .seealso: `DMMOAB`, `DMCreate()`, `DMMoabCreateVertices()`, `DMMoabCreateElement()`
1347: @*/
1348: PetscErrorCode DMMoabCreateSubmesh(DM dm, DM *newdm)
1349: {
1350: DM_Moab *dmmoab;
1351: DM_Moab *ndmmoab;
1352: moab::ErrorCode merr;
1354: PetscFunctionBegin;
1357: dmmoab = (DM_Moab *)dm->data;
1359: /* Create the basic DMMOAB object and keep the default parameters created by DM impls */
1360: PetscCall(DMMoabCreateMoab(((PetscObject)dm)->comm, dmmoab->mbiface, &dmmoab->ltog_tag, NULL, newdm));
1362: /* get all the necessary handles from the private DM object */
1363: ndmmoab = (DM_Moab *)(*newdm)->data;
1365: /* set the sub-mesh's parent DM reference */
1366: ndmmoab->parent = &dm;
1368: /* create a file set to associate all entities in current mesh */
1369: merr = ndmmoab->mbiface->create_meshset(moab::MESHSET_SET, ndmmoab->fileset);
1370: MBERR("Creating file set failed", merr);
1372: /* create a meshset and then add old fileset as child */
1373: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->vlocal);
1374: MBERR("Adding child vertices to parent failed", merr);
1375: merr = ndmmoab->mbiface->add_entities(ndmmoab->fileset, *dmmoab->elocal);
1376: MBERR("Adding child elements to parent failed", merr);
1378: /* preserve the field association between the parent and sub-mesh objects */
1379: PetscCall(DMMoabSetFieldNames(*newdm, dmmoab->numFields, dmmoab->fieldNames));
1380: PetscFunctionReturn(PETSC_SUCCESS);
1381: }
1383: PETSC_EXTERN PetscErrorCode DMMoabView_Ascii(DM dm, PetscViewer viewer)
1384: {
1385: DM_Moab *dmmoab = (DM_Moab *)(dm)->data;
1386: const char *name;
1387: MPI_Comm comm;
1388: PetscMPIInt size;
1390: PetscFunctionBegin;
1391: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1392: PetscCallMPI(MPI_Comm_size(comm, &size));
1393: PetscCall(PetscObjectGetName((PetscObject)dm, &name));
1394: PetscCall(PetscViewerASCIIPushTab(viewer));
1395: if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimensions:\n", name, dmmoab->dim));
1396: else PetscCall(PetscViewerASCIIPrintf(viewer, "Mesh in %" PetscInt_FMT " dimensions:\n", dmmoab->dim));
1397: /* print details about the global mesh */
1398: {
1399: PetscCall(PetscViewerASCIIPushTab(viewer));
1400: PetscCall(PetscViewerASCIIPrintf(viewer, "Sizes: cells=%" PetscInt_FMT ", vertices=%" PetscInt_FMT ", blocks=%" PetscInt_FMT "\n", dmmoab->nele, dmmoab->n, dmmoab->bs));
1401: /* print boundary data */
1402: PetscCall(PetscViewerASCIIPrintf(viewer, "Boundary trace:\n"));
1403: {
1404: PetscCall(PetscViewerASCIIPushTab(viewer));
1405: PetscCall(PetscViewerASCIIPrintf(viewer, "cells=%zu, faces=%zu, vertices=%zu\n", dmmoab->bndyelems->size(), dmmoab->bndyfaces->size(), dmmoab->bndyvtx->size()));
1406: PetscCall(PetscViewerASCIIPopTab(viewer));
1407: }
1408: /* print field data */
1409: PetscCall(PetscViewerASCIIPrintf(viewer, "Fields: %" PetscInt_FMT " components\n", dmmoab->numFields));
1410: {
1411: PetscCall(PetscViewerASCIIPushTab(viewer));
1412: for (int i = 0; i < dmmoab->numFields; ++i) PetscCall(PetscViewerASCIIPrintf(viewer, "[%" PetscInt_FMT "] - %s\n", i, dmmoab->fieldNames[i]));
1413: PetscCall(PetscViewerASCIIPopTab(viewer));
1414: }
1415: PetscCall(PetscViewerASCIIPopTab(viewer));
1416: }
1417: PetscCall(PetscViewerASCIIPopTab(viewer));
1418: PetscCall(PetscViewerFlush(viewer));
1419: PetscFunctionReturn(PETSC_SUCCESS);
1420: }
1422: PETSC_EXTERN PetscErrorCode DMMoabView_VTK(DM dm, PetscViewer v)
1423: {
1424: PetscFunctionReturn(PETSC_SUCCESS);
1425: }
1427: PETSC_EXTERN PetscErrorCode DMMoabView_HDF5(DM dm, PetscViewer v)
1428: {
1429: PetscFunctionReturn(PETSC_SUCCESS);
1430: }
1432: PETSC_EXTERN PetscErrorCode DMView_Moab(DM dm, PetscViewer viewer)
1433: {
1434: PetscBool iascii, ishdf5, isvtk;
1436: PetscFunctionBegin;
1439: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1440: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
1441: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
1442: if (iascii) {
1443: PetscCall(DMMoabView_Ascii(dm, viewer));
1444: } else if (ishdf5) {
1445: #if defined(PETSC_HAVE_HDF5) && defined(MOAB_HAVE_HDF5)
1446: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_VIZ));
1447: PetscCall(DMMoabView_HDF5(dm, viewer));
1448: PetscCall(PetscViewerPopFormat(viewer));
1449: #else
1450: SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
1451: #endif
1452: } else if (isvtk) {
1453: PetscCall(DMMoabView_VTK(dm, viewer));
1454: }
1455: PetscFunctionReturn(PETSC_SUCCESS);
1456: }
1458: PETSC_EXTERN PetscErrorCode DMInitialize_Moab(DM dm)
1459: {
1460: PetscFunctionBegin;
1461: dm->ops->view = DMView_Moab;
1462: dm->ops->load = NULL /* DMLoad_Moab */;
1463: dm->ops->setfromoptions = DMSetFromOptions_Moab;
1464: dm->ops->clone = DMClone_Moab;
1465: dm->ops->setup = DMSetUp_Moab;
1466: dm->ops->createlocalsection = NULL;
1467: dm->ops->createdefaultconstraints = NULL;
1468: dm->ops->createglobalvector = DMCreateGlobalVector_Moab;
1469: dm->ops->createlocalvector = DMCreateLocalVector_Moab;
1470: dm->ops->getlocaltoglobalmapping = NULL;
1471: dm->ops->createfieldis = NULL;
1472: dm->ops->createcoordinatedm = NULL /* DMCreateCoordinateDM_Moab */;
1473: dm->ops->getcoloring = NULL;
1474: dm->ops->creatematrix = DMCreateMatrix_Moab;
1475: dm->ops->createinterpolation = DMCreateInterpolation_Moab;
1476: dm->ops->createinjection = NULL /* DMCreateInjection_Moab */;
1477: dm->ops->refine = DMRefine_Moab;
1478: dm->ops->coarsen = DMCoarsen_Moab;
1479: dm->ops->refinehierarchy = DMRefineHierarchy_Moab;
1480: dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Moab;
1481: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Moab;
1482: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Moab;
1483: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Moab;
1484: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Moab;
1485: dm->ops->destroy = DMDestroy_Moab;
1486: dm->ops->createsubdm = NULL /* DMCreateSubDM_Moab */;
1487: dm->ops->getdimpoints = NULL /* DMGetDimPoints_Moab */;
1488: dm->ops->locatepoints = NULL /* DMLocatePoints_Moab */;
1489: PetscFunctionReturn(PETSC_SUCCESS);
1490: }
1492: PETSC_EXTERN PetscErrorCode DMClone_Moab(DM dm, DM *newdm)
1493: {
1494: PetscFunctionBegin;
1495: /* get all the necessary handles from the private DM object */
1496: (*newdm)->data = (DM_Moab *)dm->data;
1497: ((DM_Moab *)dm->data)->refct++;
1499: PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMMOAB));
1500: PetscCall(DMInitialize_Moab(*newdm));
1501: PetscFunctionReturn(PETSC_SUCCESS);
1502: }
1504: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)
1505: {
1506: PetscFunctionBegin;
1508: PetscCall(PetscNew((DM_Moab **)&dm->data));
1510: ((DM_Moab *)dm->data)->bs = 1;
1511: ((DM_Moab *)dm->data)->numFields = 1;
1512: ((DM_Moab *)dm->data)->n = 0;
1513: ((DM_Moab *)dm->data)->nloc = 0;
1514: ((DM_Moab *)dm->data)->nghost = 0;
1515: ((DM_Moab *)dm->data)->nele = 0;
1516: ((DM_Moab *)dm->data)->neleloc = 0;
1517: ((DM_Moab *)dm->data)->neleghost = 0;
1518: ((DM_Moab *)dm->data)->ltog_map = NULL;
1519: ((DM_Moab *)dm->data)->ltog_sendrecv = NULL;
1521: ((DM_Moab *)dm->data)->refct = 1;
1522: ((DM_Moab *)dm->data)->parent = NULL;
1523: ((DM_Moab *)dm->data)->vlocal = new moab::Range();
1524: ((DM_Moab *)dm->data)->vowned = new moab::Range();
1525: ((DM_Moab *)dm->data)->vghost = new moab::Range();
1526: ((DM_Moab *)dm->data)->elocal = new moab::Range();
1527: ((DM_Moab *)dm->data)->eghost = new moab::Range();
1529: PetscCall(DMInitialize_Moab(dm));
1530: PetscFunctionReturn(PETSC_SUCCESS);
1531: }