1: #include <petsc/private/dmmbimpl.h> /*I "petscdmmoab.h" I*/
3: #include <petscdmmoab.h>
4: #include <MBTagConventions.hpp>
5: #include <moab/Skinner.hpp>
8: /*MC
9: DMMOAB = "moab" - A DM object that encapsulates an unstructured mesh described by the MOAB mesh database.
10: Direct access to the MOAB Interface and other mesh manipulation related objects are available
11: through public API. Ability to create global and local representation of Vecs containing all
12: unknowns in the interior and shared boundary via a transparent tag-data wrapper is provided
13: along with utility functions to traverse the mesh and assemble a discrete system via
14: field-based/blocked Vec(Get/Set) methods. Input from and output to different formats are
15: available.
17: Reference: http://www.mcs.anl.gov/~fathom/moab-docs/html/contents.html
19: Level: intermediate
21: .seealso: DMType, DMMoabCreate(), DMCreate(), DMSetType(), DMMoabCreateMoab()
22: M*/
27: /*@
28: DMMoabCreate - Creates a DMMoab object, which encapsulates a moab instance
30: Collective on MPI_Comm 32: Input Parameter:
33: . comm - The communicator for the DMMoab object
35: Output Parameter:
36: . dmb - The DMMoab object
38: Level: beginner
40: .keywords: DMMoab, create
41: @*/
42: PetscErrorCodeDMMoabCreate(MPI_Comm comm, DM *dmb) 43: {
48: DMCreate(comm, dmb);
49: DMSetType(*dmb, DMMOAB);
50: return(0);
51: }
55: /*@
56: DMMoabCreate - Creates a DMMoab object, optionally from an instance and other data
58: Collective on MPI_Comm 60: Input Parameter:
61: . comm - The communicator for the DMMoab object
62: . mbiface - (ptr to) the MOAB Instance; if passed in NULL, MOAB instance is created inside PETSc, and destroyed
63: along with the DMMoab
64: . pcomm - (ptr to) a ParallelComm; if NULL, creates one internally for the whole communicator
65: . ltog_tag - A tag to use to retrieve global id for an entity; if 0, will use GLOBAL_ID_TAG_NAME/tag
66: . range - If non-NULL, contains range of entities to which DOFs will be assigned
68: Output Parameter:
69: . dmb - The DMMoab object
71: Level: intermediate
73: .keywords: DMMoab, create
74: @*/
75: PetscErrorCode DMMoabCreateMoab(MPI_Comm comm, moab::Interface *mbiface, moab::ParallelComm *pcomm, moab::Tag *ltog_tag, moab::Range *range, DM *dmb)
76: {
78: moab::ErrorCode merr;
79: moab::EntityHandle partnset;
80: PetscInt rank, nprocs;
81: DM dmmb;
82: DM_Moab *dmmoab;
87: DMMoabCreate(comm, &dmmb);
88: dmmoab = (DM_Moab*)(dmmb)->data;
90: if (!mbiface) {
91: dmmoab->mbiface = new moab::Core();
92: dmmoab->icreatedinstance = PETSC_TRUE;
93: }
94: else {
95: dmmoab->mbiface = mbiface;
96: dmmoab->icreatedinstance = PETSC_FALSE;
97: }
99: /* by default the fileset = root set. This set stores the hierarchy of entities belonging to current DM */
100: dmmoab->fileset=0;
102: if (!pcomm) {
103: MPI_Comm_rank(comm, &rank);
104: MPI_Comm_size(comm, &nprocs);
106: /* Create root sets for each mesh. Then pass these
107: to the load_file functions to be populated. */
108: merr = dmmoab->mbiface->create_meshset(moab::MESHSET_SET, partnset);MBERR("Creating partition set failed", merr);
110: /* Create the parallel communicator object with the partition handle associated with MOAB */
111: dmmoab->pcomm = moab::ParallelComm::get_pcomm(dmmoab->mbiface, partnset, &comm);
112: }
113: else {
114: DMMoabSetParallelComm(dmmb, pcomm);
115: }
117: /* do the remaining initializations for DMMoab */
118: dmmoab->bs = 1;
119: dmmoab->numFields = 1;
120: dmmoab->rw_dbglevel = 0;
121: dmmoab->partition_by_rank = PETSC_FALSE;
122: dmmoab->extra_read_options[0] = '\0';
123: dmmoab->extra_write_options[0] = '\0';
124: dmmoab->read_mode = READ_PART;
125: dmmoab->write_mode = WRITE_PART;
127: /* set global ID tag handle */
128: if (ltog_tag && *ltog_tag) {
129: DMMoabSetLocalToGlobalTag(dmmb, *ltog_tag);
130: }
131: else {
132: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);MBERRNM(merr);
133: if (ltog_tag) *ltog_tag = dmmoab->ltog_tag;
134: }
136: merr = dmmoab->mbiface->tag_get_handle(MATERIAL_SET_TAG_NAME, dmmoab->material_tag);MBERRNM(merr);
138: /* set the local range of entities (vertices) of interest */
139: if (range) {
140: DMMoabSetLocalVertices(dmmb, range);
141: }
142: *dmb=dmmb;
143: return(0);
144: }
148: /*@
149: DMMoabSetParallelComm - Set the ParallelComm used with this DMMoab
151: Collective on MPI_Comm153: Input Parameter:
154: . dm - The DMMoab object being set
155: . pcomm - The ParallelComm being set on the DMMoab
157: Level: beginner
159: .keywords: DMMoab, create
160: @*/
161: PetscErrorCodeDMMoabSetParallelComm(DM dm,moab::ParallelComm *pcomm)
162: {
163: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
168: dmmoab->pcomm = pcomm;
169: dmmoab->mbiface = pcomm->get_moab();
170: dmmoab->icreatedinstance = PETSC_FALSE;
171: return(0);
172: }
177: /*@
178: DMMoabGetParallelComm - Get the ParallelComm used with this DMMoab
180: Collective on MPI_Comm182: Input Parameter:
183: . dm - The DMMoab object being set
185: Output Parameter:
186: . pcomm - The ParallelComm for the DMMoab
188: Level: beginner
190: .keywords: DMMoab, create
191: @*/
192: PetscErrorCodeDMMoabGetParallelComm(DM dm,moab::ParallelComm **pcomm)
193: {
196: *pcomm = ((DM_Moab*)(dm)->data)->pcomm;
197: return(0);
198: }
203: /*@
204: DMMoabSetInterface - Set the MOAB instance used with this DMMoab
206: Collective on MPI_Comm208: Input Parameter:
209: . dm - The DMMoab object being set
210: . mbiface - The MOAB instance being set on this DMMoab
212: Level: beginner
214: .keywords: DMMoab, create
215: @*/
216: PetscErrorCodeDMMoabSetInterface(DM dm,moab::Interface *mbiface)
217: {
218: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
223: dmmoab->pcomm = NULL;
224: dmmoab->mbiface = mbiface;
225: dmmoab->icreatedinstance = PETSC_FALSE;
226: return(0);
227: }
232: /*@
233: DMMoabGetInterface - Get the MOAB instance used with this DMMoab
235: Collective on MPI_Comm237: Input Parameter:
238: . dm - The DMMoab object being set
240: Output Parameter:
241: . mbiface - The MOAB instance set on this DMMoab
243: Level: beginner
245: .keywords: DMMoab, create
246: @*/
247: PetscErrorCodeDMMoabGetInterface(DM dm,moab::Interface **mbiface)
248: {
249: PetscErrorCode ierr;
250: static PetscBool cite = PETSC_FALSE;
254: 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);
255: *mbiface = ((DM_Moab*)dm->data)->mbiface;
256: return(0);
257: }
262: /*@
263: DMMoabSetLocalVertices - Set the entities having DOFs on this DMMoab
265: Collective on MPI_Comm267: Input Parameter:
268: . dm - The DMMoab object being set
269: . range - The entities treated by this DMMoab
271: Level: beginner
273: .keywords: DMMoab, create
274: @*/
275: PetscErrorCodeDMMoabSetLocalVertices(DM dm,moab::Range *range)
276: {
277: moab::ErrorCode merr;
278: PetscErrorCode ierr;
279: moab::Range tmpvtxs;
280: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
284: dmmoab->vlocal->clear();
285: dmmoab->vowned->clear();
287: dmmoab->vlocal->insert(range->begin(), range->end());
289: /* filter based on parallel status */
290: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal,PSTATUS_NOT_OWNED,PSTATUS_NOT,-1,dmmoab->vowned);MBERRNM(merr);
292: /* filter all the non-owned and shared entities out of the list */
293: tmpvtxs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
294: merr = dmmoab->pcomm->filter_pstatus(tmpvtxs,PSTATUS_INTERFACE,PSTATUS_OR,-1,dmmoab->vghost);MBERRNM(merr);
295: tmpvtxs = moab::subtract(tmpvtxs, *dmmoab->vghost);
296: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, tmpvtxs);
298: /* compute and cache the sizes of local and ghosted entities */
299: dmmoab->nloc = dmmoab->vowned->size();
300: dmmoab->nghost = dmmoab->vghost->size();
301: MPI_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
302: return(0);
303: }
308: /*@
309: DMMoabGetAllVertices - Get the entities having DOFs on this DMMoab
311: Collective on MPI_Comm313: Input Parameter:
314: . dm - The DMMoab object being set
316: Output Parameter:
317: . owned - The local vertex entities in this DMMoab = (owned+ghosted)
319: Level: beginner
321: .keywords: DMMoab, create
322: @*/
323: PetscErrorCodeDMMoabGetAllVertices(DM dm,moab::Range *local)
324: {
327: if (local) *local = *((DM_Moab*)dm->data)->vlocal;
328: return(0);
329: }
335: /*@
336: DMMoabGetLocalVertices - Get the entities having DOFs on this DMMoab
338: Collective on MPI_Comm340: Input Parameter:
341: . dm - The DMMoab object being set
343: Output Parameter:
344: . owned - The owned vertex entities in this DMMoab
345: . ghost - The ghosted entities (non-owned) stored locally in this partition
347: Level: beginner
349: .keywords: DMMoab, create
350: @*/
351: PetscErrorCodeDMMoabGetLocalVertices(DM dm,const moab::Range **owned,const moab::Range **ghost)
352: {
355: if (owned) *owned = ((DM_Moab*)dm->data)->vowned;
356: if (ghost) *ghost = ((DM_Moab*)dm->data)->vghost;
357: return(0);
358: }
362: /*@
363: DMMoabGetLocalElements - Get the higher-dimensional entities that are locally owned
365: Collective on MPI_Comm367: Input Parameter:
368: . dm - The DMMoab object being set
370: Output Parameter:
371: . range - The entities owned locally
373: Level: beginner
375: .keywords: DMMoab, create
376: @*/
377: PetscErrorCodeDMMoabGetLocalElements(DM dm,const moab::Range **range)
378: {
381: if (range) *range = ((DM_Moab*)dm->data)->elocal;
382: return(0);
383: }
388: /*@
389: DMMoabSetLocalElements - Set the entities having DOFs on this DMMoab
391: Collective on MPI_Comm393: Input Parameter:
394: . dm - The DMMoab object being set
395: . range - The entities treated by this DMMoab
397: Level: beginner
399: .keywords: DMMoab, create
400: @*/
401: PetscErrorCodeDMMoabSetLocalElements(DM dm,moab::Range *range)
402: {
403: moab::ErrorCode merr;
404: PetscErrorCode ierr;
405: DM_Moab *dmmoab = (DM_Moab*)(dm)->data;
409: dmmoab->elocal->clear();
410: dmmoab->eghost->clear();
411: dmmoab->elocal->insert(range->begin(), range->end());
412: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal,PSTATUS_NOT_OWNED,PSTATUS_NOT);MBERRNM(merr);
413: *dmmoab->eghost = moab::subtract(*range, *dmmoab->elocal);
414: dmmoab->neleloc=dmmoab->elocal->size();
415: dmmoab->neleghost=dmmoab->eghost->size();
416: MPI_Allreduce(&dmmoab->nele, &dmmoab->neleloc, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
417: PetscInfo2(dm, "Created %D local and %D global elements.\n", dmmoab->neleloc, dmmoab->nele);
418: return(0);
419: }
424: /*@
425: DMMoabSetLocalToGlobalTag - Set the tag used for local to global numbering
427: Collective on MPI_Comm429: Input Parameter:
430: . dm - The DMMoab object being set
431: . ltogtag - The MOAB tag used for local to global ids
433: Level: beginner
435: .keywords: DMMoab, create
436: @*/
437: PetscErrorCodeDMMoabSetLocalToGlobalTag(DM dm,moab::Tag ltogtag)
438: {
441: ((DM_Moab*)dm->data)->ltog_tag = ltogtag;
442: return(0);
443: }
448: /*@
449: DMMoabGetLocalToGlobalTag - Get the tag used for local to global numbering
451: Collective on MPI_Comm453: Input Parameter:
454: . dm - The DMMoab object being set
456: Output Parameter:
457: . ltogtag - The MOAB tag used for local to global ids
459: Level: beginner
461: .keywords: DMMoab, create
462: @*/
463: PetscErrorCodeDMMoabGetLocalToGlobalTag(DM dm,moab::Tag *ltog_tag)
464: {
467: *ltog_tag = ((DM_Moab*)dm->data)->ltog_tag;
468: return(0);
469: }
474: /*@
475: DMMoabSetBlockSize - Set the block size used with this DMMoab
477: Collective on MPI_Comm479: Input Parameter:
480: . dm - The DMMoab object being set
481: . bs - The block size used with this DMMoab
483: Level: beginner
485: .keywords: DMMoab, create
486: @*/
487: PetscErrorCodeDMMoabSetBlockSize(DM dm,PetscInt bs)488: {
491: ((DM_Moab*)dm->data)->bs = bs;
492: return(0);
493: }
498: /*@
499: DMMoabGetBlockSize - Get the block size used with this DMMoab
501: Collective on MPI_Comm503: Input Parameter:
504: . dm - The DMMoab object being set
506: Output Parameter:
507: . bs - The block size used with this DMMoab
509: Level: beginner
511: .keywords: DMMoab, create
512: @*/
513: PetscErrorCodeDMMoabGetBlockSize(DM dm,PetscInt *bs)514: {
517: *bs = ((DM_Moab*)dm->data)->bs;
518: return(0);
519: }
524: /*@
525: DMMoabGetSize - Get the global vertex size used with this DMMoab
527: Collective on DM529: Input Parameter:
530: . dm - The DMMoab object being set
532: Output Parameter:
533: . neg - The number of global elements in the DMMoab instance
534: . nvg - The number of global vertices in the DMMoab instance
536: Level: beginner
538: .keywords: DMMoab, create
539: @*/
540: PetscErrorCodeDMMoabGetSize(DM dm,PetscInt *neg,PetscInt *nvg)541: {
544: if(neg) *neg = ((DM_Moab*)dm->data)->nele;
545: if(nvg) *nvg = ((DM_Moab*)dm->data)->n;
546: return(0);
547: }
552: /*@
553: DMMoabGetLocalSize - Get the local and ghosted vertex size used with this DMMoab
555: Collective on DM557: Input Parameter:
558: . dm - The DMMoab object being set
560: Output Parameter:
561: + nel - The number of owned elements in this processor
562: . neg - The number of ghosted elements in this processor
563: . nvl - The number of owned vertices in this processor
564: . nvg - The number of ghosted vertices in this processor
566: Level: beginner
568: .keywords: DMMoab, create
569: @*/
570: PetscErrorCodeDMMoabGetLocalSize(DM dm,PetscInt *nel,PetscInt *neg,PetscInt *nvl,PetscInt *nvg)571: {
574: if(nel) *nel = ((DM_Moab*)dm->data)->neleloc;
575: if(neg) *neg = ((DM_Moab*)dm->data)->neleghost;
576: if(nvl) *nvl = ((DM_Moab*)dm->data)->nloc;
577: if(nvg) *nvg = ((DM_Moab*)dm->data)->nghost;
578: return(0);
579: }
584: /*@
585: DMMoabGetOffset - Get the local offset for the global vector
587: Collective on MPI_Comm589: Input Parameter:
590: . dm - The DMMoab object being set
592: Output Parameter:
593: . offset - The local offset for the global vector
595: Level: beginner
597: .keywords: DMMoab, create
598: @*/
599: PetscErrorCodeDMMoabGetOffset(DM dm,PetscInt *offset)600: {
603: *offset = ((DM_Moab*)dm->data)->vstart;
604: return(0);
605: }
610: /*@
611: DMMoabGetDimension - Get the dimension of the DM Mesh
613: Collective on MPI_Comm615: Input Parameter:
616: . dm - The DMMoab object
618: Output Parameter:
619: . dim - The dimension of DM621: Level: beginner
623: .keywords: DMMoab, create
624: @*/
625: PetscErrorCodeDMMoabGetDimension(DM dm,PetscInt *dim)626: {
629: *dim = ((DM_Moab*)dm->data)->dim;
630: return(0);
631: }
636: /*@
637: DMMoabGetMaterialBlock - Get the material ID corresponding to the current entity of the DM Mesh
639: Collective on MPI_Comm641: Input Parameter:
642: . dm - The DMMoab object
643: . ehandle - The element entity handle
645: Output Parameter:
646: . mat - The material ID for the current entity
648: Level: beginner
650: .keywords: DMMoab, create
651: @*/
652: PetscErrorCodeDMMoabGetMaterialBlock(DM dm,const moab::EntityHandle ehandle, PetscInt *mat)
653: {
654: DM_Moab *dmmoab;
655: moab::ErrorCode merr;
659: if (*mat) {
660: dmmoab = (DM_Moab*)(dm)->data;
661: merr=dmmoab->mbiface->tag_get_data(dmmoab->material_tag, &ehandle, 1, mat);MBERRNM(merr);
662: }
663: return(0);
664: }
669: /*@
670: DMMoabGetVertexCoordinates - Get the coordinates corresponding to the requested vertex entities
672: Collective on MPI_Comm674: Input Parameter:
675: . dm - The DMMoab object
676: . nconn - Number of entities whose coordinates are needed
677: . conn - The vertex entity handles
679: Output Parameter:
680: . vpos - The coordinates of the requested vertex entities
682: Level: beginner
684: .seealso: DMMoabGetVertexConnectivity()
685: @*/
686: PetscErrorCodeDMMoabGetVertexCoordinates(DM dm,PetscInt nconn,const moab::EntityHandle *conn,PetscReal *vpos)
687: {
688: DM_Moab *dmmoab;
689: PetscErrorCode ierr;
690: moab::ErrorCode merr;
695: dmmoab = (DM_Moab*)(dm)->data;
697: if (!vpos) {
698: PetscMalloc1(nconn*3, &vpos);
699: }
701: /* Get connectivity information in MOAB canonical ordering */
702: merr = dmmoab->mbiface->get_coords(conn, nconn, vpos);MBERRNM(merr);
703: return(0);
704: }
709: /*@
710: DMMoabGetVertexConnectivity - Get the vertex adjacency for the given entity
712: Collective on MPI_Comm714: Input Parameter:
715: . dm - The DMMoab object
716: . vhandle - Vertex entity handle
718: Output Parameter:
719: . nconn - Number of entities whose coordinates are needed
720: . conn - The vertex entity handles
722: Level: beginner
724: .seealso: DMMoabGetVertexCoordinates(), DMMoabRestoreVertexConnectivity()
725: @*/
726: PetscErrorCodeDMMoabGetVertexConnectivity(DM dm,moab::EntityHandle vhandle,PetscInt* nconn, moab::EntityHandle **conn)
727: {
728: DM_Moab *dmmoab;
729: std::vector<moab::EntityHandle> adj_entities,connect;
730: PetscErrorCode ierr;
731: moab::ErrorCode merr;
736: dmmoab = (DM_Moab*)(dm)->data;
738: /* Get connectivity information in MOAB canonical ordering */
739: merr = dmmoab->mbiface->get_adjacencies(&vhandle, 1, 1, true, adj_entities, moab::Interface::UNION);MBERRNM(merr);
740: merr = dmmoab->mbiface->get_connectivity(&adj_entities[0],adj_entities.size(),connect);MBERRNM(merr);
742: if (conn) {
743: PetscMalloc(sizeof(moab::EntityHandle)*connect.size(), conn);
744: PetscMemcpy(*conn, &connect[0], sizeof(moab::EntityHandle)*connect.size());
745: }
746: if (nconn) *nconn=connect.size();
747: return(0);
748: }
753: /*@
754: DMMoabRestoreVertexConnectivity - Restore the vertex connectivity for the given entity
756: Collective on MPI_Comm758: Input Parameter:
759: . dm - The DMMoab object
760: . vhandle - Vertex entity handle
761: . nconn - Number of entities whose coordinates are needed
762: . conn - The vertex entity handles
764: Level: beginner
766: .seealso: DMMoabGetVertexCoordinates(), DMMoabGetVertexConnectivity()
767: @*/
768: PetscErrorCodeDMMoabRestoreVertexConnectivity(DM dm,moab::EntityHandle ehandle,PetscInt* nconn, moab::EntityHandle **conn)
769: {
770: PetscErrorCode ierr;
776: if (conn) {
777: PetscFree(*conn);
778: }
779: if (nconn) *nconn=0;
780: return(0);
781: }
786: /*@
787: DMMoabGetElementConnectivity - Get the vertex adjacency for the given entity
789: Collective on MPI_Comm791: Input Parameter:
792: . dm - The DMMoab object
793: . ehandle - Vertex entity handle
795: Output Parameter:
796: . nconn - Number of entities whose coordinates are needed
797: . conn - The vertex entity handles
799: Level: beginner
801: .seealso: DMMoabGetVertexCoordinates(), DMMoabGetVertexConnectivity(), DMMoabRestoreVertexConnectivity()
802: @*/
803: PetscErrorCodeDMMoabGetElementConnectivity(DM dm,moab::EntityHandle ehandle,PetscInt* nconn,const moab::EntityHandle **conn)
804: {
805: DM_Moab *dmmoab;
806: const moab::EntityHandle *connect;
807: moab::ErrorCode merr;
808: PetscInt nnodes;
813: dmmoab = (DM_Moab*)(dm)->data;
815: /* Get connectivity information in MOAB canonical ordering */
816: merr = dmmoab->mbiface->get_connectivity(ehandle, connect, nnodes);MBERRNM(merr);
817: if (conn) *conn=connect;
818: if (nconn) *nconn=nnodes;
819: return(0);
820: }
825: /*@
826: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
828: Collective on MPI_Comm830: Input Parameter:
831: . dm - The DMMoab object
832: . ent - Entity handle
834: Output Parameter:
835: . ent_on_boundary - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
837: Level: beginner
839: .seealso: DMMoabCheckBoundaryVertices()
840: @*/
841: PetscErrorCodeDMMoabIsEntityOnBoundary(DM dm,const moab::EntityHandle ent,PetscBool* ent_on_boundary)
842: {
843: moab::EntityType etype;
844: DM_Moab *dmmoab;
845: PetscInt edim;
850: dmmoab = (DM_Moab*)(dm)->data;
852: /* get the entity type and handle accordingly */
853: etype=dmmoab->mbiface->type_from_handle(ent);
854: if(etype >= moab::MBPOLYHEDRON) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Entity type on the boundary skin is invalid. EntityType = %D\n",etype);
856: /* get the entity dimension */
857: edim=dmmoab->mbiface->dimension_from_handle(ent);
859: *ent_on_boundary=PETSC_FALSE;
860: if(etype == moab::MBVERTEX && edim == 0) {
861: if (dmmoab->bndyvtx->index(ent) >= 0) *ent_on_boundary=PETSC_TRUE;
862: }
863: else {
864: if (edim == dmmoab->dim) { /* check the higher-dimensional elements first */
865: if (dmmoab->bndyelems->index(ent) >= 0) *ent_on_boundary=PETSC_TRUE;
866: }
867: else { /* next check the lower-dimensional faces */
868: if (dmmoab->bndyfaces->index(ent) >= 0) *ent_on_boundary=PETSC_TRUE;
869: }
870: }
871: return(0);
872: }
877: /*@
878: DMMoabIsEntityOnBoundary - Check whether a given entity is on the boundary (vertex, edge, face, element)
880: Input Parameter:
881: . dm - The DMMoab object
882: . nconn - Number of handles
883: . cnt - Array of entity handles
885: Output Parameter:
886: . isbdvtx - Array of boundary markers - PETSC_TRUE if entity on boundary; PETSC_FALSE otherwise
888: Level: beginner
890: .seealso: DMMoabIsEntityOnBoundary()
891: @*/
892: PetscErrorCode DMMoabCheckBoundaryVertices(DM dm,PetscInt nconn,const moab::EntityHandle *cnt,PetscBool* isbdvtx)
893: {
894: DM_Moab *dmmoab;
895: PetscInt i;
901: dmmoab = (DM_Moab*)(dm)->data;
903: for (i=0; i < nconn; ++i) {
904: isbdvtx[i]=(dmmoab->bndyvtx->index(cnt[i]) >= 0 ? PETSC_TRUE:PETSC_FALSE);
905: }
906: return(0);
907: }
912: /*@
913: DMMoabGetBoundaryMarkers - Return references to the vertices, faces, elements on the boundary
915: Input Parameter:
916: . dm - The DMMoab object
918: Output Parameter:
919: . bdvtx - Boundary vertices
920: . bdelems - Boundary elements
921: . bdfaces - Boundary faces
923: Level: beginner
925: .seealso: DMMoabCheckBoundaryVertices(), DMMoabIsEntityOnBoundary()
926: @*/
927: PetscErrorCodeDMMoabGetBoundaryMarkers(DM dm,const moab::Range **bdvtx,const moab::Range** bdelems,const moab::Range** bdfaces)
928: {
929: DM_Moab *dmmoab;
933: dmmoab = (DM_Moab*)(dm)->data;
935: if (bdvtx) *bdvtx = dmmoab->bndyvtx;
936: if (bdfaces) *bdfaces = dmmoab->bndyfaces;
937: if (bdelems) *bdfaces = dmmoab->bndyelems;
938: return(0);
939: }
944: PETSC_EXTERN PetscErrorCode DMDestroy_Moab(DM dm)945: {
947: PetscInt i;
948: DM_Moab *dmmoab = (DM_Moab*)dm->data;
952: if (dmmoab->icreatedinstance) {
953: delete dmmoab->mbiface;
954: }
955: dmmoab->mbiface = NULL;
956: dmmoab->pcomm = NULL;
957: delete dmmoab->vlocal;
958: delete dmmoab->vowned;
959: delete dmmoab->vghost;
960: delete dmmoab->elocal;
961: delete dmmoab->eghost;
962: delete dmmoab->bndyvtx;
963: delete dmmoab->bndyfaces;
964: delete dmmoab->bndyelems;
966: PetscFree(dmmoab->gsindices);
967: PetscFree2(dmmoab->gidmap,dmmoab->lidmap);
968: PetscFree2(dmmoab->lgmap,dmmoab->llmap);
969: PetscFree(dmmoab->dfill);
970: PetscFree(dmmoab->ofill);
971: if (dmmoab->fieldNames) {
972: for(i=0; i<dmmoab->numFields; i++) {
973: PetscFree(dmmoab->fieldNames[i]);
974: }
975: PetscFree(dmmoab->fieldNames);
976: }
977: VecScatterDestroy(&dmmoab->ltog_sendrecv);
978: ISLocalToGlobalMappingDestroy(&dmmoab->ltog_map);
979: PetscFree(dm->data);
980: return(0);
981: }
986: PETSC_EXTERN PetscErrorCode DMSetFromOptions_Moab(PetscOptions *PetscOptionsObject,DM dm)987: {
989: DM_Moab *dmmoab = (DM_Moab*)dm->data;
993: PetscOptionsHead(PetscOptionsObject,"DMMoab Options");
994: PetscOptionsInt("-dm_moab_rw_dbg", "The verbosity level for reading and writing MOAB meshes", "DMView", dmmoab->rw_dbglevel, &dmmoab->rw_dbglevel, NULL);
995: 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);
996: /* TODO: typically, the read options are needed before a DM is completely created and available in which case, the options wont be available ?? */
997: PetscOptionsString("-dm_moab_read_opts", "Extra options to enable MOAB reader to load DM from file", "DMView", dmmoab->extra_read_options, dmmoab->extra_read_options, PETSC_MAX_PATH_LEN, NULL);
998: PetscOptionsString("-dm_moab_write_opts", "Extra options to enable MOAB writer to serialize DM to file", "DMView", dmmoab->extra_write_options, dmmoab->extra_write_options, PETSC_MAX_PATH_LEN, NULL);
999: PetscOptionsEnum("-dm_moab_read_mode", "MOAB parallel read mode", "DMView", MoabReadModes, (PetscEnum)dmmoab->read_mode, (PetscEnum*)&dmmoab->read_mode, NULL);
1000: PetscOptionsEnum("-dm_moab_write_mode", "MOAB parallel write mode", "DMView", MoabWriteModes, (PetscEnum)dmmoab->write_mode, (PetscEnum*)&dmmoab->write_mode, NULL);
1001: return(0);
1002: }
1007: PETSC_EXTERN PetscErrorCode DMSetUp_Moab(DM dm)1008: {
1009: PetscErrorCode ierr;
1010: moab::ErrorCode merr;
1011: Vec local, global;
1012: IS from,to;
1013: moab::Range::iterator iter;
1014: PetscInt i,j,f,bs,gmin,lmin,lmax,vent,totsize;
1015: DM_Moab *dmmoab = (DM_Moab*)dm->data;
1016: moab::Range adjs;
1020: /* Get the local and shared vertices and cache it */
1021: if (dmmoab->mbiface == PETSC_NULL || dmmoab->pcomm == PETSC_NULL) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ORDER, "Set the MOAB Interface and ParallelComm objects before calling SetUp.");
1023: /* Get the entities recursively in the current part of the mesh, if user did not set the local vertices explicitly */
1024: if (dmmoab->vlocal->empty())
1025: {
1026: merr = dmmoab->mbiface->get_entities_by_type(dmmoab->fileset,moab::MBVERTEX,*dmmoab->vlocal,true);MBERRNM(merr);
1028: /* filter based on parallel status */
1029: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->vlocal,PSTATUS_NOT_OWNED,PSTATUS_NOT,-1,dmmoab->vowned);MBERRNM(merr);
1031: /* filter all the non-owned and shared entities out of the list */
1032: adjs = moab::subtract(*dmmoab->vlocal, *dmmoab->vowned);
1033: merr = dmmoab->pcomm->filter_pstatus(adjs,PSTATUS_INTERFACE,PSTATUS_OR,-1,dmmoab->vghost);MBERRNM(merr);
1034: adjs = moab::subtract(adjs, *dmmoab->vghost);
1035: *dmmoab->vlocal = moab::subtract(*dmmoab->vlocal, adjs);
1037: /* compute and cache the sizes of local and ghosted entities */
1038: dmmoab->nloc = dmmoab->vowned->size();
1039: dmmoab->nghost = dmmoab->vghost->size();
1040: MPI_Allreduce(&dmmoab->nloc, &dmmoab->n, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
1041: }
1043: {
1044: /* get the information about the local elements in the mesh */
1045: dmmoab->eghost->clear();
1047: /* first decipher the leading dimension */
1048: for (i=3;i>0;i--) {
1049: dmmoab->elocal->clear();
1050: merr = dmmoab->mbiface->get_entities_by_dimension(dmmoab->fileset, i, *dmmoab->elocal, true);CHKERRQ(merr);
1052: /* store the current mesh dimension */
1053: if (dmmoab->elocal->size()) {
1054: dmmoab->dim=i;
1055: break;
1056: }
1057: }
1059: /* filter the ghosted and owned element list */
1060: *dmmoab->eghost = *dmmoab->elocal;
1061: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->elocal,PSTATUS_NOT_OWNED,PSTATUS_NOT);MBERRNM(merr);
1062: *dmmoab->eghost = moab::subtract(*dmmoab->eghost, *dmmoab->elocal);
1064: dmmoab->neleloc = dmmoab->elocal->size();
1065: dmmoab->neleghost = dmmoab->eghost->size();
1066: MPI_Allreduce(&dmmoab->neleloc, &dmmoab->nele, 1, MPI_INTEGER, MPI_SUM, ((PetscObject)dm)->comm);
1067: }
1069: bs = dmmoab->bs;
1070: if (!dmmoab->ltog_tag) {
1071: /* Get the global ID tag. The global ID tag is applied to each
1072: vertex. It acts as an global identifier which MOAB uses to
1073: assemble the individual pieces of the mesh */
1074: merr = dmmoab->mbiface->tag_get_handle(GLOBAL_ID_TAG_NAME, dmmoab->ltog_tag);MBERRNM(merr);
1075: }
1077: totsize=dmmoab->vlocal->size();
1078: if (totsize != dmmoab->nloc+dmmoab->nghost) SETERRQ2(PETSC_COMM_WORLD, PETSC_ERR_ARG_OUTOFRANGE, "Mismatch between local and owned+ghost vertices. %D != %D.",totsize,dmmoab->nloc+dmmoab->nghost);
1079: PetscMalloc1(totsize,&dmmoab->gsindices);
1080: {
1081: /* first get the local indices */
1082: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag,*dmmoab->vowned,&dmmoab->gsindices[0]);MBERRNM(merr);
1083: /* next get the ghosted indices */
1084: if (dmmoab->nghost) {
1085: merr = dmmoab->mbiface->tag_get_data(dmmoab->ltog_tag,*dmmoab->vghost,&dmmoab->gsindices[dmmoab->nloc]);MBERRNM(merr);
1086: }
1088: /* find out the local and global minima of GLOBAL_ID */
1089: lmin=lmax=dmmoab->gsindices[0];
1090: for (i=0; i<totsize; ++i) {
1091: if(lmin>dmmoab->gsindices[i]) lmin=dmmoab->gsindices[i];
1092: if(lmax<dmmoab->gsindices[i]) lmax=dmmoab->gsindices[i];
1093: }
1095: MPI_Allreduce(&lmin, &gmin, 1, MPI_INT, MPI_MIN, ((PetscObject)dm)->comm);
1097: /* set the GID map */
1098: for (i=0; i<totsize; ++i) {
1099: dmmoab->gsindices[i]-=gmin; /* zero based index needed for IS */
1100: }
1101: lmin-=gmin;
1102: lmax-=gmin;
1104: PetscInfo3(NULL, "GLOBAL_ID: Local minima - %D, Local maxima - %D, Global minima - %D.\n", lmin, lmax, gmin);
1105: }
1106: if(!(dmmoab->bs==dmmoab->numFields || dmmoab->bs == 1)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mismatch between block size and number of component fields. %D != 1 OR %D != %D.",dmmoab->bs,dmmoab->bs,dmmoab->numFields);
1108: {
1109: i=(PetscInt)dmmoab->vlocal->back()+1;
1110: //i=(PetscInt)(dmmoab->vlocal->back()-dmmoab->vlocal->front())+1;
1111: j=totsize*dmmoab->numFields;
1112: PetscMalloc2(i,&dmmoab->gidmap,i,&dmmoab->lidmap);
1113: PetscMalloc2(j,&dmmoab->lgmap,j,&dmmoab->llmap);
1115: i=j=0;
1116: /* set the owned vertex data first */
1117: for(moab::Range::iterator iter = dmmoab->vowned->begin(); iter != dmmoab->vowned->end(); iter++,i++) {
1118: vent=(PetscInt)(*iter);
1119: dmmoab->gidmap[vent]=dmmoab->gsindices[i];
1120: dmmoab->lidmap[vent]=i;
1121: if (bs > 1) {
1122: for (f=0;f<dmmoab->numFields;f++,j++) {
1123: dmmoab->lgmap[j]=dmmoab->gsindices[i]*dmmoab->numFields+f;
1124: dmmoab->llmap[j]=i*dmmoab->numFields+f;
1125: }
1126: }
1127: else {
1128: for (f=0;f<dmmoab->numFields;f++,j++) {
1129: dmmoab->lgmap[j]=totsize*f+dmmoab->gsindices[i];
1130: dmmoab->llmap[j]=totsize*f+i;
1131: }
1132: }
1133: }
1134: /* next arrange all the ghosted data information */
1135: for(moab::Range::iterator iter = dmmoab->vghost->begin(); iter != dmmoab->vghost->end(); iter++,i++) {
1136: vent=(PetscInt)(*iter);
1137: dmmoab->gidmap[vent]=dmmoab->gsindices[i];
1138: dmmoab->lidmap[vent]=i;
1139: if (bs > 1) {
1140: for (f=0;f<dmmoab->numFields;f++,j++) {
1141: dmmoab->lgmap[j]=dmmoab->gsindices[i]*dmmoab->numFields+f;
1142: dmmoab->llmap[j]=i*dmmoab->numFields+f;
1143: }
1144: }
1145: else {
1146: for (f=0;f<dmmoab->numFields;f++,j++) {
1147: dmmoab->lgmap[j]=totsize*f+dmmoab->gsindices[i];
1148: dmmoab->llmap[j]=totsize*f+i;
1149: }
1150: }
1151: }
1153: /* We need to create the Global to Local Vector Scatter Contexts
1154: 1) First create a local and global vector
1155: 2) Create a local and global IS1156: 3) Create VecScatter and LtoGMapping objects
1157: 4) Cleanup the IS and Vec objects
1158: */
1159: DMCreateGlobalVector(dm, &global);
1160: DMCreateLocalVector(dm, &local);
1162: VecGetOwnershipRange(global, &dmmoab->vstart, &dmmoab->vend);
1163: PetscInfo3(NULL, "Total-size = %D\t Owned = %D, Ghosted = %D.\n", totsize, dmmoab->nloc, dmmoab->nghost);
1165: /* global to local must retrieve ghost points */
1166: ISCreateStride(((PetscObject)dm)->comm,dmmoab->nloc*dmmoab->numFields,dmmoab->vstart,1,&from);
1167: ISSetBlockSize(from,bs);
1169: ISCreateGeneral(((PetscObject)dm)->comm,dmmoab->nloc*dmmoab->numFields,&dmmoab->lgmap[0],PETSC_COPY_VALUES,&to);
1170: ISSetBlockSize(to,bs);
1172: if (!dmmoab->ltog_map) {
1173: /* create to the local to global mapping for vectors in order to use VecSetValuesLocal */
1174: ISLocalToGlobalMappingCreate(((PetscObject)dm)->comm,dmmoab->bs,totsize*dmmoab->numFields,dmmoab->lgmap,
1175: PETSC_COPY_VALUES,&dmmoab->ltog_map);
1176: }
1178: /* now create the scatter object from local to global vector */
1179: VecScatterCreate(local,from,global,to,&dmmoab->ltog_sendrecv);
1181: /* clean up IS, Vec */
1182: ISDestroy(&from);
1183: ISDestroy(&to);
1184: VecDestroy(&local);
1185: VecDestroy(&global);
1186: }
1188: /* skin the boundary and store nodes */
1189: {
1190: /* get the skin vertices of boundary faces for the current partition and then filter
1191: the local, boundary faces, vertices and elements alone via PSTATUS flags;
1192: this should not give us any ghosted boundary, but if user needs such a functionality
1193: it would be easy to add it based on the find_skin query below */
1194: moab::Skinner skinner(dmmoab->mbiface);
1196: dmmoab->bndyvtx = new moab::Range();
1197: dmmoab->bndyfaces = new moab::Range();
1198: dmmoab->bndyelems = new moab::Range();
1200: /* get the entities on the skin - only the faces */
1201: merr = skinner.find_skin(dmmoab->fileset, *dmmoab->elocal, false, *dmmoab->bndyfaces, NULL, false, true, false);MBERRNM(merr); // 'false' param indicates we want faces back, not vertices
1203: /* filter all the non-owned and shared entities out of the list */
1204: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces,PSTATUS_NOT_OWNED,PSTATUS_NOT);MBERRNM(merr);
1205: merr = dmmoab->pcomm->filter_pstatus(*dmmoab->bndyfaces,PSTATUS_SHARED,PSTATUS_NOT);MBERRNM(merr);
1207: /* get all the nodes via connectivity and the parent elements via adjacency information */
1208: merr = dmmoab->mbiface->get_connectivity(*dmmoab->bndyfaces, *dmmoab->bndyvtx, false);MBERRNM(ierr);
1209: merr = dmmoab->mbiface->get_adjacencies(*dmmoab->bndyfaces, dmmoab->dim, false, *dmmoab->bndyelems, moab::Interface::UNION);MBERRNM(ierr);
1210: PetscInfo3(NULL, "Found %D boundary vertices, %D boundary faces and %D boundary elements.\n", dmmoab->bndyvtx->size(), dmmoab->bndyvtx->size(), dmmoab->bndyelems->size());
1211: }
1212: return(0);
1213: }
1217: PETSC_EXTERN PetscErrorCode DMCreate_Moab(DM dm)1218: {
1223: PetscNewLog(dm,(DM_Moab**)&dm->data);
1225: ((DM_Moab*)dm->data)->bs = 1;
1226: ((DM_Moab*)dm->data)->numFields = 1;
1227: ((DM_Moab*)dm->data)->n = 0;
1228: ((DM_Moab*)dm->data)->nloc = 0;
1229: ((DM_Moab*)dm->data)->nghost = 0;
1230: ((DM_Moab*)dm->data)->nele = 0;
1231: ((DM_Moab*)dm->data)->neleloc = 0;
1232: ((DM_Moab*)dm->data)->neleghost = 0;
1233: ((DM_Moab*)dm->data)->ltog_map = PETSC_NULL;
1234: ((DM_Moab*)dm->data)->ltog_sendrecv = PETSC_NULL;
1236: ((DM_Moab*)dm->data)->vlocal = new moab::Range();
1237: ((DM_Moab*)dm->data)->vowned = new moab::Range();
1238: ((DM_Moab*)dm->data)->vghost = new moab::Range();
1239: ((DM_Moab*)dm->data)->elocal = new moab::Range();
1240: ((DM_Moab*)dm->data)->eghost = new moab::Range();
1242: dm->ops->createglobalvector = DMCreateGlobalVector_Moab;
1243: dm->ops->createlocalvector = DMCreateLocalVector_Moab;
1244: dm->ops->creatematrix = DMCreateMatrix_Moab;
1245: dm->ops->setup = DMSetUp_Moab;
1246: dm->ops->destroy = DMDestroy_Moab;
1247: dm->ops->setfromoptions = DMSetFromOptions_Moab;
1248: dm->ops->globaltolocalbegin = DMGlobalToLocalBegin_Moab;
1249: dm->ops->globaltolocalend = DMGlobalToLocalEnd_Moab;
1250: dm->ops->localtoglobalbegin = DMLocalToGlobalBegin_Moab;
1251: dm->ops->localtoglobalend = DMLocalToGlobalEnd_Moab;
1252: return(0);
1253: }