Actual source code: dmmbvec.cxx

      1: #include <petsc/private/dmmbimpl.h>
      2: #include <petsc/private/vecimpl.h>

      4: #include <petscdmmoab.h>
      5: #include <MBTagConventions.hpp>

      7: #define USE_NATIVE_PETSCVEC

      9: /* declare some private DMMoab specific overrides */
     10: static PetscErrorCode DMCreateVector_Moab_Private(DM dm,moab::Tag tag,const moab::Range* userrange,PetscBool is_global_vec,PetscBool destroy_tag,Vec *vec);
     11: static PetscErrorCode DMVecUserDestroy_Moab(void *user);
     12: static PetscErrorCode DMVecDuplicate_Moab(Vec x,Vec *y);
     13: #ifdef MOAB_HAVE_MPI
     14: static PetscErrorCode DMVecCreateTagName_Moab_Private(moab::Interface *mbiface,moab::ParallelComm *pcomm,char** tag_name);
     15: #else
     16: static PetscErrorCode DMVecCreateTagName_Moab_Private(moab::Interface *mbiface,char** tag_name);
     17: #endif

     19: /*@C
     20:   DMMoabCreateVector - Create a Vec from either an existing tag, or a specified tag size, and a range of entities

     22:   Collective

     24:   Input Parameter:
     25: + dm              - The DMMoab object being set
     26: . tag             - If non-zero, block size will be taken from the tag size
     27: . range           - If non-empty, Vec corresponds to these entities, otherwise to the entities set on the DMMoab
     28: . is_global_vec   - If true, this is a local representation of the Vec (including ghosts in parallel), otherwise a truly parallel one
     29: - destroy_tag     - If true, MOAB tag is destroyed with Vec, otherwise it is left on MOAB

     31:   Output Parameter:
     32: . vec             - The created vector

     34:   Level: beginner

     36: .seealso: VecCreate()
     37: @*/
     38: PetscErrorCode DMMoabCreateVector(DM dm, moab::Tag tag, const moab::Range* range, PetscBool is_global_vec, PetscBool destroy_tag, Vec *vec)
     39: {
     40:   PetscErrorCode     ierr;

     43:   if (!tag && (!range || range->empty())) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Both tag and range cannot be null.");

     45:   DMCreateVector_Moab_Private(dm, tag, range, is_global_vec, destroy_tag, vec);
     46:   return(0);
     47: }

     50: /*@C
     51:   DMMoabGetVecTag - Get the MOAB tag associated with this Vec

     53:   Input Parameter:
     54: . vec - Vec being queried

     56:   Output Parameter:
     57: . tag - Tag associated with this Vec. NULL if vec is a native PETSc Vec.

     59:   Level: beginner

     61: .seealso: DMMoabCreateVector(), DMMoabGetVecRange()
     62: @*/
     63: PetscErrorCode DMMoabGetVecTag(Vec vec, moab::Tag *tag)
     64: {
     65:   PetscContainer  moabdata;
     66:   Vec_MOAB        *vmoab;
     67:   PetscErrorCode  ierr;

     72:   /* Get the MOAB private data */
     73:   PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject*) &moabdata);
     74:   PetscContainerGetPointer(moabdata, (void**) &vmoab);

     76:   *tag = vmoab->tag;
     77:   return(0);
     78: }

     81: /*@C
     82:   DMMoabGetVecRange - Get the MOAB entities associated with this Vec

     84:   Input Parameter:
     85: . vec   - Vec being queried

     87:   Output Parameter:
     88: . range - Entities associated with this Vec. NULL if vec is a native PETSc Vec.

     90:   Level: beginner

     92: .seealso: DMMoabCreateVector(), DMMoabGetVecTag()
     93: @*/
     94: PetscErrorCode DMMoabGetVecRange(Vec vec, moab::Range *range)
     95: {
     96:   PetscContainer  moabdata;
     97:   Vec_MOAB        *vmoab;
     98:   PetscErrorCode  ierr;

    103:   /* Get the MOAB private data handle */
    104:   PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject*) &moabdata);
    105:   PetscContainerGetPointer(moabdata, (void**) &vmoab);

    107:   *range = *vmoab->tag_range;
    108:   return(0);
    109: }

    112: /*@C
    113:   DMMoabVecGetArray - Returns the writable direct access array to the local representation of MOAB tag data for the underlying vector using locally owned+ghosted range of entities

    115:   Collective

    117:   Input Parameter:
    118: + dm              - The DMMoab object being set
    119: - vec             - The Vector whose underlying data is requested

    121:   Output Parameter:
    122: . array           - The local data array

    124:   Level: intermediate

    126: .seealso: DMMoabVecRestoreArray(), DMMoabVecGetArrayRead(), DMMoabVecRestoreArrayRead()
    127: @*/
    128: PetscErrorCode  DMMoabVecGetArray(DM dm, Vec vec, void* array)
    129: {
    130:   DM_Moab        *dmmoab;
    131:   moab::ErrorCode merr;
    132:   PetscErrorCode  ierr;
    133:   PetscInt        count, i, f;
    134:   moab::Tag       vtag;
    135:   PetscScalar     **varray;
    136:   PetscScalar     *marray;
    137:   PetscContainer  moabdata;
    138:   Vec_MOAB        *vmoab, *xmoab;

    144:   dmmoab = (DM_Moab*)dm->data;

    146:   /* Get the Vec_MOAB struct for the original vector */
    147:   PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject*) &moabdata);
    148:   PetscContainerGetPointer(moabdata, (void**)&vmoab);

    150:   /* Get the real scalar array handle */
    151:   varray = reinterpret_cast<PetscScalar**>(array);

    153:   if (vmoab->is_native_vec) {

    155:     /* get the local representation of the arrays from Vectors */
    156:     DMGetLocalVector(dm, &vmoab->local);
    157:     DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vmoab->local);
    158:     DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vmoab->local);

    160:     /* Get the Vec_MOAB struct for the original vector */
    161:     PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject*) &moabdata);
    162:     PetscContainerGetPointer(moabdata, (void**)&xmoab);

    164:     /* get the local representation of the arrays from Vectors */
    165:     VecGhostGetLocalForm(vmoab->local, &xmoab->local);
    166:     VecGhostUpdateBegin(vmoab->local, INSERT_VALUES, SCATTER_FORWARD);
    167:     VecGhostUpdateEnd(vmoab->local, INSERT_VALUES, SCATTER_FORWARD);

    169:     VecGetArray(xmoab->local, varray);
    170:   }
    171:   else {

    173:     /* Get the MOAB private data */
    174:     DMMoabGetVecTag(vec, &vtag);

    176: #ifdef MOAB_HAVE_MPI
    177:     /* exchange the data into ghost cells first */
    178:     merr = dmmoab->pcomm->exchange_tags(vtag, *dmmoab->vlocal); MBERRNM(merr);
    179: #endif

    181:     PetscMalloc1((dmmoab->nloc + dmmoab->nghost) * dmmoab->numFields, varray);

    183:     /* Get the array data for local entities */
    184:     merr = dmmoab->mbiface->tag_iterate(vtag, dmmoab->vlocal->begin(), dmmoab->vlocal->end(), count, reinterpret_cast<void*&>(marray), false); MBERRNM(merr);
    185:     if (count != dmmoab->nloc + dmmoab->nghost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mismatch between local vertices and tag partition for Vec. %D != %D.", count, dmmoab->nloc + dmmoab->nghost);

    187:     i = 0;
    188:     for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
    189:       for (f = 0; f < dmmoab->numFields; f++, i++)
    190:         (*varray)[dmmoab->lidmap[(PetscInt)*iter - dmmoab->seqstart]*dmmoab->numFields + f] = marray[i];
    191:       //(*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]]=marray[i];
    192:     }
    193:   }
    194:   return(0);
    195: }

    198: /*@C
    199:   DMMoabVecRestoreArray - Restores the writable direct access array obtained via DMMoabVecGetArray

    201:   Collective

    203:   Input Parameter:
    204: + dm              - The DMMoab object being set
    205: . vec             - The Vector whose underlying data is requested
    206: - array           - The local data array

    208:   Level: intermediate

    210: .seealso: DMMoabVecGetArray(), DMMoabVecGetArrayRead(), DMMoabVecRestoreArrayRead()
    211: @*/
    212: PetscErrorCode  DMMoabVecRestoreArray(DM dm, Vec vec, void* array)
    213: {
    214:   DM_Moab        *dmmoab;
    215:   moab::ErrorCode merr;
    216:   PetscErrorCode  ierr;
    217:   moab::Tag       vtag;
    218:   PetscInt        count, i, f;
    219:   PetscScalar     **varray;
    220:   PetscScalar     *marray;
    221:   PetscContainer  moabdata;
    222:   Vec_MOAB        *vmoab, *xmoab;

    228:   dmmoab = (DM_Moab*)dm->data;

    230:   /* Get the Vec_MOAB struct for the original vector */
    231:   PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject*) &moabdata);
    232:   PetscContainerGetPointer(moabdata, (void**)&vmoab);

    234:   /* Get the real scalar array handle */
    235:   varray = reinterpret_cast<PetscScalar**>(array);

    237:   if (vmoab->is_native_vec) {

    239:     /* Get the Vec_MOAB struct for the original vector */
    240:     PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject*) &moabdata);
    241:     PetscContainerGetPointer(moabdata, (void**)&xmoab);

    243:     /* get the local representation of the arrays from Vectors */
    244:     VecRestoreArray(xmoab->local, varray);
    245:     VecGhostUpdateBegin(vmoab->local, ADD_VALUES, SCATTER_REVERSE);
    246:     VecGhostUpdateEnd(vmoab->local, ADD_VALUES, SCATTER_REVERSE);
    247:     VecGhostRestoreLocalForm(vmoab->local, &xmoab->local);

    249:     /* restore local pieces */
    250:     DMLocalToGlobalBegin(dm, vmoab->local, INSERT_VALUES, vec);
    251:     DMLocalToGlobalEnd(dm, vmoab->local, INSERT_VALUES, vec);
    252:     DMRestoreLocalVector(dm, &vmoab->local);
    253:   }
    254:   else {

    256:     /* Get the MOAB private data */
    257:     DMMoabGetVecTag(vec, &vtag);

    259:     /* Get the array data for local entities */
    260:     merr = dmmoab->mbiface->tag_iterate(vtag, dmmoab->vlocal->begin(), dmmoab->vlocal->end(), count, reinterpret_cast<void*&>(marray), false); MBERRNM(merr);
    261:     if (count != dmmoab->nloc + dmmoab->nghost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mismatch between local vertices and tag partition for Vec. %D != %D.", count, dmmoab->nloc + dmmoab->nghost);

    263:     i = 0;
    264:     for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
    265:       for (f = 0; f < dmmoab->numFields; f++, i++)
    266:         marray[i] = (*varray)[dmmoab->lidmap[(PetscInt) * iter - dmmoab->seqstart] * dmmoab->numFields + f];
    267:       //marray[i] = (*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]];
    268:     }

    270: #ifdef MOAB_HAVE_MPI
    271:     /* reduce the tags correctly -> should probably let the user choose how to reduce in the future
    272:       For all FEM residual based assembly calculations, MPI_SUM should serve well */
    273:     merr = dmmoab->pcomm->reduce_tags(vtag, MPI_SUM, *dmmoab->vlocal); MBERRV(dmmoab->mbiface, merr);
    274: #endif
    275:     PetscFree(*varray);
    276:   }
    277:   return(0);
    278: }

    280: /*@C
    281:   DMMoabVecGetArrayRead - Returns the read-only direct access array to the local representation of MOAB tag data for the underlying vector using locally owned+ghosted range of entities

    283:   Collective

    285:   Input Parameter:
    286: + dm              - The DMMoab object being set
    287: - vec             - The Vector whose underlying data is requested

    289:   Output Parameter:
    290: . array           - The local data array

    292:   Level: intermediate

    294: .seealso: DMMoabVecRestoreArrayRead(), DMMoabVecGetArray(), DMMoabVecRestoreArray()
    295: @*/
    296: PetscErrorCode  DMMoabVecGetArrayRead(DM dm, Vec vec, void* array)
    297: {
    298:   DM_Moab        *dmmoab;
    299:   moab::ErrorCode merr;
    300:   PetscErrorCode  ierr;
    301:   PetscInt        count, i, f;
    302:   moab::Tag       vtag;
    303:   PetscScalar     **varray;
    304:   PetscScalar     *marray;
    305:   PetscContainer  moabdata;
    306:   Vec_MOAB        *vmoab, *xmoab;

    312:   dmmoab = (DM_Moab*)dm->data;

    314:   /* Get the Vec_MOAB struct for the original vector */
    315:   PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject*) &moabdata);
    316:   PetscContainerGetPointer(moabdata, (void**)&vmoab);

    318:   /* Get the real scalar array handle */
    319:   varray = reinterpret_cast<PetscScalar**>(array);

    321:   if (vmoab->is_native_vec) {
    322:     /* get the local representation of the arrays from Vectors */
    323:     DMGetLocalVector(dm, &vmoab->local);
    324:     DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vmoab->local);
    325:     DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vmoab->local);

    327:     /* Get the Vec_MOAB struct for the original vector */
    328:     PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject*) &moabdata);
    329:     PetscContainerGetPointer(moabdata, (void**)&xmoab);

    331:     /* get the local representation of the arrays from Vectors */
    332:     VecGhostGetLocalForm(vmoab->local, &xmoab->local);
    333:     VecGhostUpdateBegin(vmoab->local, INSERT_VALUES, SCATTER_FORWARD);
    334:     VecGhostUpdateEnd(vmoab->local, INSERT_VALUES, SCATTER_FORWARD);
    335:     VecGetArray(xmoab->local, varray);
    336:   }
    337:   else {
    338:     /* Get the MOAB private data */
    339:     DMMoabGetVecTag(vec, &vtag);

    341: #ifdef MOAB_HAVE_MPI
    342:     /* exchange the data into ghost cells first */
    343:     merr = dmmoab->pcomm->exchange_tags(vtag, *dmmoab->vlocal); MBERRNM(merr);
    344: #endif
    345:     PetscMalloc1((dmmoab->nloc + dmmoab->nghost) * dmmoab->numFields, varray);

    347:     /* Get the array data for local entities */
    348:     merr = dmmoab->mbiface->tag_iterate(vtag, dmmoab->vlocal->begin(), dmmoab->vlocal->end(), count, reinterpret_cast<void*&>(marray), false); MBERRNM(merr);
    349:     if (count != dmmoab->nloc + dmmoab->nghost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Mismatch between local vertices and tag partition for Vec. %D != %D.", count, dmmoab->nloc + dmmoab->nghost);

    351:     i = 0;
    352:     for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
    353:       for (f = 0; f < dmmoab->numFields; f++, i++)
    354:         (*varray)[dmmoab->lidmap[(PetscInt)*iter - dmmoab->seqstart]*dmmoab->numFields + f] = marray[i];
    355:       //(*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]]=marray[i];
    356:     }
    357:   }
    358:   return(0);
    359: }

    362: /*@C
    363:   DMMoabVecRestoreArray - Restores the read-only direct access array obtained via DMMoabVecGetArray

    365:   Collective

    367:   Input Parameter:
    368: + dm              - The DMMoab object being set
    369: . vec             - The Vector whose underlying data is requested
    370: - array           - The local data array

    372:   Level: intermediate

    374: .seealso: DMMoabVecGetArrayRead(), DMMoabVecGetArray(), DMMoabVecRestoreArray()
    375: @*/
    376: PetscErrorCode  DMMoabVecRestoreArrayRead(DM dm, Vec vec, void* array)
    377: {
    378:   PetscErrorCode  ierr;
    379:   PetscScalar     **varray;
    380:   PetscContainer  moabdata;
    381:   Vec_MOAB        *vmoab, *xmoab;

    388:   /* Get the Vec_MOAB struct for the original vector */
    389:   PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject*) &moabdata);
    390:   PetscContainerGetPointer(moabdata, (void**)&vmoab);

    392:   /* Get the real scalar array handle */
    393:   varray = reinterpret_cast<PetscScalar**>(array);

    395:   if (vmoab->is_native_vec) {
    396:     /* Get the Vec_MOAB struct for the original vector */
    397:     PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject*) &moabdata);
    398:     PetscContainerGetPointer(moabdata, (void**)&xmoab);

    400:     /* restore the local representation of the arrays from Vectors */
    401:     VecRestoreArray(xmoab->local, varray);
    402:     VecGhostRestoreLocalForm(vmoab->local, &xmoab->local);

    404:     /* restore local pieces */
    405:     DMRestoreLocalVector(dm, &vmoab->local);
    406:   }
    407:   else {
    408:     /* Nothing to do but just free the memory allocated before */
    409:     PetscFree(*varray);

    411:   }
    412:   return(0);
    413: }

    416: PetscErrorCode DMCreateVector_Moab_Private(DM dm, moab::Tag tag, const moab::Range* userrange, PetscBool is_global_vec, PetscBool destroy_tag, Vec *vec)
    417: {
    418:   PetscErrorCode    ierr;
    419:   moab::ErrorCode   merr;
    420:   PetscBool         is_newtag;
    421:   const moab::Range *range;
    422:   PetscInt          count, lnative_vec, gnative_vec;
    423:   std::string       ttname;
    424:   PetscScalar       *data_ptr, *defaultvals;

    426:   Vec_MOAB *vmoab;
    427:   DM_Moab *dmmoab = (DM_Moab*)dm->data;
    428: #ifdef MOAB_HAVE_MPI
    429:   moab::ParallelComm *pcomm = dmmoab->pcomm;
    430: #endif
    431:   moab::Interface *mbiface = dmmoab->mbiface;

    434:   if (sizeof(PetscReal) != sizeof(PetscScalar)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "MOAB tags only support Real types (Complex-type unsupported)");
    435:   if (!userrange) range = dmmoab->vowned;
    436:   else range = userrange;
    437:   if (!range) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Input range cannot be empty or call DMSetUp first.");

    439: #ifndef USE_NATIVE_PETSCVEC
    440:   /* If the tag data is in a single sequence, use PETSc native vector since tag_iterate isn't useful anymore */
    441:   lnative_vec = (range->psize() - 1);
    442: #else
    443:   lnative_vec = 1; /* NOTE: Testing PETSc vector: will force to create native vector all the time */
    444: //  lnative_vec=0; /* NOTE: Testing MOAB vector: will force to create MOAB tag_iterate based vector all the time */
    445: #endif

    447: #ifdef MOAB_HAVE_MPI
    448:   MPIU_Allreduce(&lnative_vec, &gnative_vec, 1, MPI_INT, MPI_MAX, (((PetscObject)dm)->comm));
    449: #else
    450:   gnative_vec = lnative_vec;
    451: #endif

    453:   /* Create the MOAB internal data object */
    454:   PetscNew(&vmoab);
    455:   vmoab->is_native_vec = (gnative_vec > 0 ? PETSC_TRUE : PETSC_FALSE);

    457:   if (!vmoab->is_native_vec) {
    458:     merr = moab::MB_SUCCESS;
    459:     if (tag != 0) merr = mbiface->tag_get_name(tag, ttname);
    460:     if (!ttname.length() || merr != moab::MB_SUCCESS) {
    461:       /* get the new name for the anonymous MOABVec -> the tag_name will be destroyed along with Tag */
    462:       char *tag_name = NULL;
    463: #ifdef MOAB_HAVE_MPI
    464:       DMVecCreateTagName_Moab_Private(mbiface,pcomm,&tag_name);
    465: #else
    466:       DMVecCreateTagName_Moab_Private(mbiface,&tag_name);
    467: #endif
    468:       is_newtag = PETSC_TRUE;

    470:       /* Create the default value for the tag (all zeros) */
    471:       PetscCalloc1(dmmoab->numFields, &defaultvals);

    473:       /* Create the tag */
    474:       merr = mbiface->tag_get_handle(tag_name, dmmoab->numFields, moab::MB_TYPE_DOUBLE, tag,
    475:                                      moab::MB_TAG_DENSE | moab::MB_TAG_CREAT, defaultvals); MBERRNM(merr);
    476:       PetscFree(tag_name);
    477:       PetscFree(defaultvals);
    478:     }
    479:     else {
    480:       /* Make sure the tag data is of type "double" */
    481:       moab::DataType tag_type;
    482:       merr = mbiface->tag_get_data_type(tag, tag_type); MBERRNM(merr);
    483:       if (tag_type != moab::MB_TYPE_DOUBLE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Tag data type must be MB_TYPE_DOUBLE");
    484:       is_newtag = destroy_tag;
    485:     }

    487:     vmoab->tag = tag;
    488:     vmoab->new_tag = is_newtag;
    489:   }
    490:   vmoab->mbiface = mbiface;
    491: #ifdef MOAB_HAVE_MPI
    492:   vmoab->pcomm = pcomm;
    493: #endif
    494:   vmoab->is_global_vec = is_global_vec;
    495:   vmoab->tag_size = dmmoab->bs;

    497:   if (vmoab->is_native_vec) {

    499:     /* Create the PETSc Vector directly and attach our functions accordingly */
    500:     if (!is_global_vec) {
    501:       /* This is an MPI Vector with ghosted padding */
    502:       VecCreateGhostBlock((((PetscObject)dm)->comm), dmmoab->bs, dmmoab->numFields * dmmoab->nloc,
    503:                                  dmmoab->numFields * dmmoab->n, dmmoab->nghost, &dmmoab->gsindices[dmmoab->nloc], vec);
    504:     }
    505:     else {
    506:       /* This is an MPI/SEQ Vector */
    507:       VecCreate((((PetscObject)dm)->comm), vec);
    508:       VecSetSizes(*vec, dmmoab->numFields * dmmoab->nloc, PETSC_DECIDE);
    509:       VecSetBlockSize(*vec, dmmoab->bs);
    510:       VecSetType(*vec, VECMPI);
    511:     }
    512:   }
    513:   else {
    514:     /* Call tag_iterate. This will cause MOAB to allocate memory for the
    515:        tag data if it hasn't already happened */
    516:     merr = mbiface->tag_iterate(tag, range->begin(), range->end(), count, (void*&)data_ptr); MBERRNM(merr);

    518:     /* set the reference for vector range */
    519:     vmoab->tag_range = new moab::Range(*range);
    520:     merr = mbiface->tag_get_length(tag, dmmoab->numFields); MBERRNM(merr);

    522:     /* Create the PETSc Vector
    523:       Query MOAB mesh to check if there are any ghosted entities
    524:         -> if we do, create a ghosted vector to map correctly to the same layout
    525:         -> else, create a non-ghosted parallel vector */
    526:     if (!is_global_vec) {
    527:       /* This is an MPI Vector with ghosted padding */
    528:       VecCreateGhostBlockWithArray((((PetscObject)dm)->comm), dmmoab->bs, dmmoab->numFields * dmmoab->nloc,
    529:                                           dmmoab->numFields * dmmoab->n, dmmoab->nghost, &dmmoab->gsindices[dmmoab->nloc], data_ptr, vec);
    530:     }
    531:     else {
    532:       /* This is an MPI Vector without ghosted padding */
    533:       VecCreateMPIWithArray((((PetscObject)dm)->comm), dmmoab->bs, dmmoab->numFields * range->size(),
    534:                                    PETSC_DECIDE, data_ptr, vec);
    535:     }
    536:   }
    537:   VecSetFromOptions(*vec);

    539:   /* create a container and store the internal MOAB data for faster access based on Entities etc */
    540:   PetscContainer moabdata;
    541:   PetscContainerCreate(PETSC_COMM_WORLD, &moabdata);
    542:   PetscContainerSetPointer(moabdata, vmoab);
    543:   PetscContainerSetUserDestroy(moabdata, DMVecUserDestroy_Moab);
    544:   PetscObjectCompose((PetscObject) * vec, "MOABData", (PetscObject)moabdata);
    545:   (*vec)->ops->duplicate = DMVecDuplicate_Moab;
    546:   PetscContainerDestroy(&moabdata);

    548:   /* Vector created, manually set local to global mapping */
    549:   if (dmmoab->ltog_map) {
    550:     VecSetLocalToGlobalMapping(*vec, dmmoab->ltog_map);
    551:   }

    553:   /* set the DM reference to the vector */
    554:   VecSetDM(*vec, dm);
    555:   return(0);
    556: }

    559: /*  DMVecCreateTagName_Moab_Private
    560:  *
    561:  *  Creates a unique tag name that will be shared across processes. If
    562:  *  pcomm is NULL, then this is a serial vector. A unique tag name
    563:  *  will be returned in tag_name in either case.
    564:  *
    565:  *  The tag names have the format _PETSC_VEC_N where N is some integer.
    566:  *
    567:  *  NOTE: The tag_name is allocated in this routine; The user needs to free
    568:  *        the character array.
    569:  */
    570: #ifdef MOAB_HAVE_MPI
    571: PetscErrorCode DMVecCreateTagName_Moab_Private(moab::Interface *mbiface, moab::ParallelComm *pcomm, char** tag_name)
    572: #else
    573: PetscErrorCode DMVecCreateTagName_Moab_Private(moab::Interface *mbiface, char** tag_name)
    574: #endif
    575: {
    576:   moab::ErrorCode mberr;
    577:   PetscErrorCode  ierr;
    578:   PetscInt        n, global_n;
    579:   moab::Tag indexTag;

    582:   const char*       PVEC_PREFIX      = "__PETSC_VEC_";
    583:   PetscMalloc1(PETSC_MAX_PATH_LEN, tag_name);

    585:   moab::EntityHandle rootset = mbiface->get_root_set();

    587:   /* Check to see if there are any PETSc vectors defined */
    588:   /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
    589:   mberr = mbiface->tag_get_handle("__PETSC_VECS__", 1, moab::MB_TYPE_INTEGER, indexTag,
    590:                                   moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT, 0); MBERRNM(mberr);
    591:   mberr = mbiface->tag_get_data(indexTag, &rootset, 1, &n);
    592:   if (mberr == moab::MB_TAG_NOT_FOUND) n = 0; /* this is the first temporary vector */
    593:   else MBERRNM(mberr);

    595:   /* increment the new value of n */
    596:   ++n;

    598: #ifdef MOAB_HAVE_MPI
    599:   /* Make sure that n is consistent across all processes */
    600:   MPIU_Allreduce(&n, &global_n, 1, MPI_INT, MPI_MAX, pcomm->comm());
    601: #else
    602:   global_n = n;
    603: #endif

    605:   /* Set the new name accordingly and return */
    606:   PetscSNPrintf(*tag_name, PETSC_MAX_PATH_LEN - 1, "%s_%D", PVEC_PREFIX, global_n);
    607:   mberr = mbiface->tag_set_data(indexTag, &rootset, 1, (const void*)&global_n); MBERRNM(mberr);
    608:   return(0);
    609: }

    612: PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM dm, Vec *gvec)
    613: {
    614:   PetscErrorCode  ierr;
    615:   DM_Moab         *dmmoab = (DM_Moab*)dm->data;

    620:   DMCreateVector_Moab_Private(dm,NULL,dmmoab->vowned,PETSC_TRUE,PETSC_TRUE,gvec);
    621:   return(0);
    622: }

    625: PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM dm, Vec *lvec)
    626: {
    627:   PetscErrorCode  ierr;
    628:   DM_Moab         *dmmoab = (DM_Moab*)dm->data;

    633:   DMCreateVector_Moab_Private(dm,NULL,dmmoab->vlocal,PETSC_FALSE,PETSC_TRUE,lvec);
    634:   return(0);
    635: }

    638: PetscErrorCode DMVecDuplicate_Moab(Vec x, Vec *y)
    639: {
    641:   DM             dm;
    642:   PetscContainer  moabdata;
    643:   Vec_MOAB        *vmoab;

    649:   /* Get the Vec_MOAB struct for the original vector */
    650:   PetscObjectQuery((PetscObject)x, "MOABData", (PetscObject*) &moabdata);
    651:   PetscContainerGetPointer(moabdata, (void**)&vmoab);

    653:   VecGetDM(x, &dm);

    656:   DMCreateVector_Moab_Private(dm,NULL,vmoab->tag_range,vmoab->is_global_vec,PETSC_TRUE,y);
    657:   VecSetDM(*y, dm);
    658:   return(0);
    659: }

    662: PetscErrorCode DMVecUserDestroy_Moab(void *user)
    663: {
    664:   Vec_MOAB        *vmoab = (Vec_MOAB*)user;
    665:   PetscErrorCode  ierr;
    666:   moab::ErrorCode merr;

    669:   if (vmoab->new_tag && vmoab->tag) {
    670:     /* Tag was created via a call to VecDuplicate, delete the underlying tag in MOAB */
    671:     merr = vmoab->mbiface->tag_delete(vmoab->tag); MBERRNM(merr);
    672:   }
    673:   delete vmoab->tag_range;
    674:   vmoab->tag = NULL;
    675:   vmoab->mbiface = NULL;
    676: #ifdef MOAB_HAVE_MPI
    677:   vmoab->pcomm = NULL;
    678: #endif
    679:   PetscFree(vmoab);
    680:   return(0);
    681: }

    684: PETSC_EXTERN PetscErrorCode  DMGlobalToLocalBegin_Moab(DM dm, Vec g, InsertMode mode, Vec l)
    685: {
    686:   PetscErrorCode    ierr;
    687:   DM_Moab         *dmmoab = (DM_Moab*)dm->data;

    690:   VecScatterBegin(dmmoab->ltog_sendrecv, g, l, mode, SCATTER_REVERSE);
    691:   return(0);
    692: }

    695: PETSC_EXTERN PetscErrorCode  DMGlobalToLocalEnd_Moab(DM dm, Vec g, InsertMode mode, Vec l)
    696: {
    697:   PetscErrorCode    ierr;
    698:   DM_Moab         *dmmoab = (DM_Moab*)dm->data;

    701:   VecScatterEnd(dmmoab->ltog_sendrecv, g, l, mode, SCATTER_REVERSE);
    702:   return(0);
    703: }

    706: PETSC_EXTERN PetscErrorCode  DMLocalToGlobalBegin_Moab(DM dm, Vec l, InsertMode mode, Vec g)
    707: {
    708:   PetscErrorCode    ierr;
    709:   DM_Moab         *dmmoab = (DM_Moab*)dm->data;

    712:   VecScatterBegin(dmmoab->ltog_sendrecv, l, g, mode, SCATTER_FORWARD);
    713:   return(0);
    714: }

    717: PETSC_EXTERN PetscErrorCode  DMLocalToGlobalEnd_Moab(DM dm, Vec l, InsertMode mode, Vec g)
    718: {
    719:   PetscErrorCode    ierr;
    720:   DM_Moab         *dmmoab = (DM_Moab*)dm->data;

    723:   VecScatterEnd(dmmoab->ltog_sendrecv, l, g, mode, SCATTER_FORWARD);
    724:   return(0);
    725: }