Actual source code: dmmoab.cxx

petsc-3.6.1 2015-08-06
Report Typos and Errors
  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: PetscErrorCode DMMoabCreate(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_Comm

153:   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: PetscErrorCode DMMoabSetParallelComm(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_Comm

182:   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: PetscErrorCode DMMoabGetParallelComm(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_Comm

208:   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: PetscErrorCode DMMoabSetInterface(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_Comm

237:   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: PetscErrorCode DMMoabGetInterface(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_Comm

267:   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: PetscErrorCode DMMoabSetLocalVertices(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_Comm

313:   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: PetscErrorCode DMMoabGetAllVertices(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_Comm

340:   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: PetscErrorCode DMMoabGetLocalVertices(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_Comm

367:   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: PetscErrorCode DMMoabGetLocalElements(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_Comm

393:   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: PetscErrorCode DMMoabSetLocalElements(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_Comm

429:   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: PetscErrorCode DMMoabSetLocalToGlobalTag(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_Comm

453:   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: PetscErrorCode DMMoabGetLocalToGlobalTag(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_Comm

479:   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: PetscErrorCode DMMoabSetBlockSize(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_Comm

503:   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: PetscErrorCode DMMoabGetBlockSize(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 DM

529:   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: PetscErrorCode DMMoabGetSize(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 DM

557:   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: PetscErrorCode DMMoabGetLocalSize(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_Comm

589:   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: PetscErrorCode DMMoabGetOffset(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_Comm

615:   Input Parameter:
616: . dm - The DMMoab object

618:   Output Parameter:
619: . dim - The dimension of DM

621:   Level: beginner

623: .keywords: DMMoab, create
624: @*/
625: PetscErrorCode DMMoabGetDimension(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_Comm

641:   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: PetscErrorCode DMMoabGetMaterialBlock(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_Comm

674:   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: PetscErrorCode DMMoabGetVertexCoordinates(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_Comm

714:   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: PetscErrorCode DMMoabGetVertexConnectivity(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_Comm

758:   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: PetscErrorCode DMMoabRestoreVertexConnectivity(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_Comm

791:   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: PetscErrorCode DMMoabGetElementConnectivity(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_Comm

830:   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: PetscErrorCode DMMoabIsEntityOnBoundary(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: PetscErrorCode DMMoabGetBoundaryMarkers(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 IS
1156:        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: }