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 Parameters:
     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:   PetscFunctionBegin;
     41:   PetscCheck(tag || (range && !range->empty()), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Both tag and range cannot be null.");

     43:   PetscCall(DMCreateVector_Moab_Private(dm, tag, range, is_global_vec, destroy_tag, vec));
     44:   PetscFunctionReturn(PETSC_SUCCESS);
     45: }

     47: /*@C
     48:   DMMoabGetVecTag - Get the MOAB tag associated with this Vec

     50:   Input Parameter:
     51: . vec - Vec being queried

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

     56:   Level: beginner

     58: .seealso: `DMMoabCreateVector()`, `DMMoabGetVecRange()`
     59: @*/
     60: PetscErrorCode DMMoabGetVecTag(Vec vec, moab::Tag *tag)
     61: {
     62:   PetscContainer moabdata;
     63:   Vec_MOAB      *vmoab;

     65:   PetscFunctionBegin;
     66:   PetscAssertPointer(tag, 2);

     68:   /* Get the MOAB private data */
     69:   PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
     70:   PetscCall(PetscContainerGetPointer(moabdata, (void **)&vmoab));

     72:   *tag = vmoab->tag;
     73:   PetscFunctionReturn(PETSC_SUCCESS);
     74: }

     76: /*@C
     77:   DMMoabGetVecRange - Get the MOAB entities associated with this Vec

     79:   Input Parameter:
     80: . vec - Vec being queried

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

     85:   Level: beginner

     87: .seealso: `DMMoabCreateVector()`, `DMMoabGetVecTag()`
     88: @*/
     89: PetscErrorCode DMMoabGetVecRange(Vec vec, moab::Range *range)
     90: {
     91:   PetscContainer moabdata;
     92:   Vec_MOAB      *vmoab;

     94:   PetscFunctionBegin;
     95:   PetscAssertPointer(range, 2);

     97:   /* Get the MOAB private data handle */
     98:   PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
     99:   PetscCall(PetscContainerGetPointer(moabdata, (void **)&vmoab));

    101:   *range = *vmoab->tag_range;
    102:   PetscFunctionReturn(PETSC_SUCCESS);
    103: }

    105: /*@C
    106:   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

    108:   Collective

    110:   Input Parameters:
    111: + dm  - The DMMoab object being set
    112: - vec - The Vector whose underlying data is requested

    114:   Output Parameter:
    115: . array - The local data array

    117:   Level: intermediate

    119: .seealso: `DMMoabVecRestoreArray()`, `DMMoabVecGetArrayRead()`, `DMMoabVecRestoreArrayRead()`
    120: @*/
    121: PetscErrorCode DMMoabVecGetArray(DM dm, Vec vec, void *array)
    122: {
    123:   DM_Moab        *dmmoab;
    124:   moab::ErrorCode merr;
    125:   PetscInt        count, i, f;
    126:   moab::Tag       vtag;
    127:   PetscScalar   **varray;
    128:   PetscScalar    *marray;
    129:   PetscContainer  moabdata;
    130:   Vec_MOAB       *vmoab, *xmoab;

    132:   PetscFunctionBegin;
    135:   PetscAssertPointer(array, 3);
    136:   dmmoab = (DM_Moab *)dm->data;

    138:   /* Get the Vec_MOAB struct for the original vector */
    139:   PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
    140:   PetscCall(PetscContainerGetPointer(moabdata, (void **)&vmoab));

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

    145:   if (vmoab->is_native_vec) {
    146:     /* get the local representation of the arrays from Vectors */
    147:     PetscCall(DMGetLocalVector(dm, &vmoab->local));
    148:     PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vmoab->local));
    149:     PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vmoab->local));

    151:     /* Get the Vec_MOAB struct for the original vector */
    152:     PetscCall(PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject *)&moabdata));
    153:     PetscCall(PetscContainerGetPointer(moabdata, (void **)&xmoab));

    155:     /* get the local representation of the arrays from Vectors */
    156:     PetscCall(VecGhostGetLocalForm(vmoab->local, &xmoab->local));
    157:     PetscCall(VecGhostUpdateBegin(vmoab->local, INSERT_VALUES, SCATTER_FORWARD));
    158:     PetscCall(VecGhostUpdateEnd(vmoab->local, INSERT_VALUES, SCATTER_FORWARD));

    160:     PetscCall(VecGetArray(xmoab->local, varray));
    161:   } else {
    162:     /* Get the MOAB private data */
    163:     PetscCall(DMMoabGetVecTag(vec, &vtag));

    165: #ifdef MOAB_HAVE_MPI
    166:     /* exchange the data into ghost cells first */
    167:     merr = dmmoab->pcomm->exchange_tags(vtag, *dmmoab->vlocal);
    168:     MBERRNM(merr);
    169: #endif

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

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

    178:     i = 0;
    179:     for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
    180:       for (f = 0; f < dmmoab->numFields; f++, i++) (*varray)[dmmoab->lidmap[(PetscInt)*iter - dmmoab->seqstart] * dmmoab->numFields + f] = marray[i];
    181:       //(*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]]=marray[i];
    182:     }
    183:   }
    184:   PetscFunctionReturn(PETSC_SUCCESS);
    185: }

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

    190:   Collective

    192:   Input Parameters:
    193: + dm    - The DMMoab object being set
    194: . vec   - The Vector whose underlying data is requested
    195: - array - The local data array

    197:   Level: intermediate

    199: .seealso: `DMMoabVecGetArray()`, `DMMoabVecGetArrayRead()`, `DMMoabVecRestoreArrayRead()`
    200: @*/
    201: PetscErrorCode DMMoabVecRestoreArray(DM dm, Vec vec, void *array)
    202: {
    203:   DM_Moab        *dmmoab;
    204:   moab::ErrorCode merr;
    205:   moab::Tag       vtag;
    206:   PetscInt        count, i, f;
    207:   PetscScalar   **varray;
    208:   PetscScalar    *marray;
    209:   PetscContainer  moabdata;
    210:   Vec_MOAB       *vmoab, *xmoab;

    212:   PetscFunctionBegin;
    215:   PetscAssertPointer(array, 3);
    216:   dmmoab = (DM_Moab *)dm->data;

    218:   /* Get the Vec_MOAB struct for the original vector */
    219:   PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
    220:   PetscCall(PetscContainerGetPointer(moabdata, (void **)&vmoab));

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

    225:   if (vmoab->is_native_vec) {
    226:     /* Get the Vec_MOAB struct for the original vector */
    227:     PetscCall(PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject *)&moabdata));
    228:     PetscCall(PetscContainerGetPointer(moabdata, (void **)&xmoab));

    230:     /* get the local representation of the arrays from Vectors */
    231:     PetscCall(VecRestoreArray(xmoab->local, varray));
    232:     PetscCall(VecGhostUpdateBegin(vmoab->local, ADD_VALUES, SCATTER_REVERSE));
    233:     PetscCall(VecGhostUpdateEnd(vmoab->local, ADD_VALUES, SCATTER_REVERSE));
    234:     PetscCall(VecGhostRestoreLocalForm(vmoab->local, &xmoab->local));

    236:     /* restore local pieces */
    237:     PetscCall(DMLocalToGlobalBegin(dm, vmoab->local, INSERT_VALUES, vec));
    238:     PetscCall(DMLocalToGlobalEnd(dm, vmoab->local, INSERT_VALUES, vec));
    239:     PetscCall(DMRestoreLocalVector(dm, &vmoab->local));
    240:   } else {
    241:     /* Get the MOAB private data */
    242:     PetscCall(DMMoabGetVecTag(vec, &vtag));

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

    249:     i = 0;
    250:     for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
    251:       for (f = 0; f < dmmoab->numFields; f++, i++) marray[i] = (*varray)[dmmoab->lidmap[(PetscInt)*iter - dmmoab->seqstart] * dmmoab->numFields + f];
    252:       //marray[i] = (*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]];
    253:     }

    255: #ifdef MOAB_HAVE_MPI
    256:     /* reduce the tags correctly -> should probably let the user choose how to reduce in the future
    257:       For all FEM residual based assembly calculations, MPI_SUM should serve well */
    258:     merr = dmmoab->pcomm->reduce_tags(vtag, MPI_SUM, *dmmoab->vlocal);
    259:     MBERRV(dmmoab->mbiface, merr);
    260: #endif
    261:     PetscCall(PetscFree(*varray));
    262:   }
    263:   PetscFunctionReturn(PETSC_SUCCESS);
    264: }

    266: /*@C
    267:   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

    269:   Collective

    271:   Input Parameters:
    272: + dm  - The DMMoab object being set
    273: - vec - The Vector whose underlying data is requested

    275:   Output Parameter:
    276: . array - The local data array

    278:   Level: intermediate

    280: .seealso: `DMMoabVecRestoreArrayRead()`, `DMMoabVecGetArray()`, `DMMoabVecRestoreArray()`
    281: @*/
    282: PetscErrorCode DMMoabVecGetArrayRead(DM dm, Vec vec, void *array)
    283: {
    284:   DM_Moab        *dmmoab;
    285:   moab::ErrorCode merr;
    286:   PetscInt        count, i, f;
    287:   moab::Tag       vtag;
    288:   PetscScalar   **varray;
    289:   PetscScalar    *marray;
    290:   PetscContainer  moabdata;
    291:   Vec_MOAB       *vmoab, *xmoab;

    293:   PetscFunctionBegin;
    296:   PetscAssertPointer(array, 3);
    297:   dmmoab = (DM_Moab *)dm->data;

    299:   /* Get the Vec_MOAB struct for the original vector */
    300:   PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
    301:   PetscCall(PetscContainerGetPointer(moabdata, (void **)&vmoab));

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

    306:   if (vmoab->is_native_vec) {
    307:     /* get the local representation of the arrays from Vectors */
    308:     PetscCall(DMGetLocalVector(dm, &vmoab->local));
    309:     PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vmoab->local));
    310:     PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vmoab->local));

    312:     /* Get the Vec_MOAB struct for the original vector */
    313:     PetscCall(PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject *)&moabdata));
    314:     PetscCall(PetscContainerGetPointer(moabdata, (void **)&xmoab));

    316:     /* get the local representation of the arrays from Vectors */
    317:     PetscCall(VecGhostGetLocalForm(vmoab->local, &xmoab->local));
    318:     PetscCall(VecGhostUpdateBegin(vmoab->local, INSERT_VALUES, SCATTER_FORWARD));
    319:     PetscCall(VecGhostUpdateEnd(vmoab->local, INSERT_VALUES, SCATTER_FORWARD));
    320:     PetscCall(VecGetArray(xmoab->local, varray));
    321:   } else {
    322:     /* Get the MOAB private data */
    323:     PetscCall(DMMoabGetVecTag(vec, &vtag));

    325: #ifdef MOAB_HAVE_MPI
    326:     /* exchange the data into ghost cells first */
    327:     merr = dmmoab->pcomm->exchange_tags(vtag, *dmmoab->vlocal);
    328:     MBERRNM(merr);
    329: #endif
    330:     PetscCall(PetscMalloc1((dmmoab->nloc + dmmoab->nghost) * dmmoab->numFields, varray));

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

    337:     i = 0;
    338:     for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
    339:       for (f = 0; f < dmmoab->numFields; f++, i++) (*varray)[dmmoab->lidmap[(PetscInt)*iter - dmmoab->seqstart] * dmmoab->numFields + f] = marray[i];
    340:       //(*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]]=marray[i];
    341:     }
    342:   }
    343:   PetscFunctionReturn(PETSC_SUCCESS);
    344: }

    346: /*@C
    347:   DMMoabVecRestoreArrayRead - Restores the read-only direct access array obtained via DMMoabVecGetArray

    349:   Collective

    351:   Input Parameters:
    352: + dm    - The DMMoab object being set
    353: . vec   - The Vector whose underlying data is requested
    354: - array - The local data array

    356:   Level: intermediate

    358: .seealso: `DMMoabVecGetArrayRead()`, `DMMoabVecGetArray()`, `DMMoabVecRestoreArray()`
    359: @*/
    360: PetscErrorCode DMMoabVecRestoreArrayRead(DM dm, Vec vec, void *array)
    361: {
    362:   PetscScalar  **varray;
    363:   PetscContainer moabdata;
    364:   Vec_MOAB      *vmoab, *xmoab;

    366:   PetscFunctionBegin;
    369:   PetscAssertPointer(array, 3);

    371:   /* Get the Vec_MOAB struct for the original vector */
    372:   PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
    373:   PetscCall(PetscContainerGetPointer(moabdata, (void **)&vmoab));

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

    378:   if (vmoab->is_native_vec) {
    379:     /* Get the Vec_MOAB struct for the original vector */
    380:     PetscCall(PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject *)&moabdata));
    381:     PetscCall(PetscContainerGetPointer(moabdata, (void **)&xmoab));

    383:     /* restore the local representation of the arrays from Vectors */
    384:     PetscCall(VecRestoreArray(xmoab->local, varray));
    385:     PetscCall(VecGhostRestoreLocalForm(vmoab->local, &xmoab->local));

    387:     /* restore local pieces */
    388:     PetscCall(DMRestoreLocalVector(dm, &vmoab->local));
    389:   } else {
    390:     /* Nothing to do but just free the memory allocated before */
    391:     PetscCall(PetscFree(*varray));
    392:   }
    393:   PetscFunctionReturn(PETSC_SUCCESS);
    394: }

    396: static PetscErrorCode DMCreateVector_Moab_Private(DM dm, moab::Tag tag, const moab::Range *userrange, PetscBool is_global_vec, PetscBool destroy_tag, Vec *vec)
    397: {
    398:   moab::ErrorCode    merr;
    399:   PetscBool          is_newtag;
    400:   const moab::Range *range;
    401:   PetscInt           count, lnative_vec, gnative_vec;
    402:   std::string        ttname;
    403:   PetscScalar       *data_ptr, *defaultvals;

    405:   Vec_MOAB *vmoab;
    406:   DM_Moab  *dmmoab = (DM_Moab *)dm->data;
    407: #ifdef MOAB_HAVE_MPI
    408:   moab::ParallelComm *pcomm = dmmoab->pcomm;
    409: #endif
    410:   moab::Interface *mbiface = dmmoab->mbiface;

    412:   PetscFunctionBegin;
    413:   PetscCheck(sizeof(PetscReal) == sizeof(PetscScalar), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "MOAB tags only support Real types (Complex-type unsupported)");
    414:   if (!userrange) range = dmmoab->vowned;
    415:   else range = userrange;
    416:   PetscCheck(range, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Input range cannot be empty or call DMSetUp first.");

    418: #ifndef USE_NATIVE_PETSCVEC
    419:   /* If the tag data is in a single sequence, use PETSc native vector since tag_iterate isn't useful anymore */
    420:   lnative_vec = (range->psize() - 1);
    421: #else
    422:   lnative_vec = 1; /* NOTE: Testing PETSc vector: will force to create native vector all the time */
    423:                    //  lnative_vec=0; /* NOTE: Testing MOAB vector: will force to create MOAB tag_iterate based vector all the time */
    424: #endif

    426: #ifdef MOAB_HAVE_MPI
    427:   PetscCallMPI(MPIU_Allreduce(&lnative_vec, &gnative_vec, 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm));
    428: #else
    429:   gnative_vec = lnative_vec;
    430: #endif

    432:   /* Create the MOAB internal data object */
    433:   PetscCall(PetscNew(&vmoab));
    434:   vmoab->is_native_vec = (gnative_vec > 0 ? PETSC_TRUE : PETSC_FALSE);

    436:   if (!vmoab->is_native_vec) {
    437:     merr = moab::MB_SUCCESS;
    438:     if (tag != 0) merr = mbiface->tag_get_name(tag, ttname);
    439:     if (!ttname.length() || merr != moab::MB_SUCCESS) {
    440:       /* get the new name for the anonymous MOABVec -> the tag_name will be destroyed along with Tag */
    441:       char *tag_name = NULL;
    442: #ifdef MOAB_HAVE_MPI
    443:       PetscCall(DMVecCreateTagName_Moab_Private(mbiface, pcomm, &tag_name));
    444: #else
    445:       PetscCall(DMVecCreateTagName_Moab_Private(mbiface, &tag_name));
    446: #endif
    447:       is_newtag = PETSC_TRUE;

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

    452:       /* Create the tag */
    453:       merr = mbiface->tag_get_handle(tag_name, dmmoab->numFields, moab::MB_TYPE_DOUBLE, tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT, defaultvals);
    454:       MBERRNM(merr);
    455:       PetscCall(PetscFree(tag_name));
    456:       PetscCall(PetscFree(defaultvals));
    457:     } else {
    458:       /* Make sure the tag data is of type "double" */
    459:       moab::DataType tag_type;
    460:       merr = mbiface->tag_get_data_type(tag, tag_type);
    461:       MBERRNM(merr);
    462:       PetscCheck(tag_type == moab::MB_TYPE_DOUBLE, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Tag data type must be MB_TYPE_DOUBLE");
    463:       is_newtag = destroy_tag;
    464:     }

    466:     vmoab->tag     = tag;
    467:     vmoab->new_tag = is_newtag;
    468:   }
    469:   vmoab->mbiface = mbiface;
    470: #ifdef MOAB_HAVE_MPI
    471:   vmoab->pcomm = pcomm;
    472: #endif
    473:   vmoab->is_global_vec = is_global_vec;
    474:   vmoab->tag_size      = dmmoab->bs;

    476:   if (vmoab->is_native_vec) {
    477:     /* Create the PETSc Vector directly and attach our functions accordingly */
    478:     if (!is_global_vec) {
    479:       /* This is an MPI Vector with ghosted padding */
    480:       PetscCall(VecCreateGhostBlock(((PetscObject)dm)->comm, dmmoab->bs, dmmoab->numFields * dmmoab->nloc, dmmoab->numFields * dmmoab->n, dmmoab->nghost, &dmmoab->gsindices[dmmoab->nloc], vec));
    481:     } else {
    482:       /* This is an MPI/SEQ Vector */
    483:       PetscCall(VecCreate(((PetscObject)dm)->comm, vec));
    484:       PetscCall(VecSetSizes(*vec, dmmoab->numFields * dmmoab->nloc, PETSC_DECIDE));
    485:       PetscCall(VecSetBlockSize(*vec, dmmoab->bs));
    486:       PetscCall(VecSetType(*vec, VECMPI));
    487:     }
    488:   } else {
    489:     /* Call tag_iterate. This will cause MOAB to allocate memory for the
    490:        tag data if it hasn't already happened */
    491:     merr = mbiface->tag_iterate(tag, range->begin(), range->end(), count, (void *&)data_ptr);
    492:     MBERRNM(merr);

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

    499:     /* Create the PETSc Vector
    500:       Query MOAB mesh to check if there are any ghosted entities
    501:         -> if we do, create a ghosted vector to map correctly to the same layout
    502:         -> else, create a non-ghosted parallel vector */
    503:     if (!is_global_vec) {
    504:       /* This is an MPI Vector with ghosted padding */
    505:       PetscCall(VecCreateGhostBlockWithArray(((PetscObject)dm)->comm, dmmoab->bs, dmmoab->numFields * dmmoab->nloc, dmmoab->numFields * dmmoab->n, dmmoab->nghost, &dmmoab->gsindices[dmmoab->nloc], data_ptr, vec));
    506:     } else {
    507:       /* This is an MPI Vector without ghosted padding */
    508:       PetscCall(VecCreateMPIWithArray(((PetscObject)dm)->comm, dmmoab->bs, dmmoab->numFields * range->size(), PETSC_DECIDE, data_ptr, vec));
    509:     }
    510:   }
    511:   PetscCall(VecSetFromOptions(*vec));

    513:   /* create a container and store the internal MOAB data for faster access based on Entities etc */
    514:   PetscCall(PetscObjectContainerCompose((PetscObject)*vec, "MOABData", vmoab, DMVecUserDestroy_Moab));
    515:   (*vec)->ops->duplicate = DMVecDuplicate_Moab;

    517:   /* Vector created, manually set local to global mapping */
    518:   if (dmmoab->ltog_map) PetscCall(VecSetLocalToGlobalMapping(*vec, dmmoab->ltog_map));

    520:   /* set the DM reference to the vector */
    521:   PetscCall(VecSetDM(*vec, dm));
    522:   PetscFunctionReturn(PETSC_SUCCESS);
    523: }

    525: /*  DMVecCreateTagName_Moab_Private
    526:  *
    527:  *  Creates a unique tag name that will be shared across processes. If
    528:  *  pcomm is NULL, then this is a serial vector. A unique tag name
    529:  *  will be returned in tag_name in either case.
    530:  *
    531:  *  The tag names have the format _PETSC_VEC_N where N is some integer.
    532:  *
    533:  *  NOTE: The tag_name is allocated in this routine; The user needs to free
    534:  *        the character array.
    535:  */
    536: #ifdef MOAB_HAVE_MPI
    537: PetscErrorCode DMVecCreateTagName_Moab_Private(moab::Interface *mbiface, moab::ParallelComm *pcomm, char **tag_name)
    538: #else
    539: static PetscErrorCode DMVecCreateTagName_Moab_Private(moab::Interface *mbiface, char **tag_name)
    540: #endif
    541: {
    542:   moab::ErrorCode mberr;
    543:   PetscInt        n, global_n;
    544:   moab::Tag       indexTag;

    546:   PetscFunctionBegin;
    547:   const char *PVEC_PREFIX = "__PETSC_VEC_";
    548:   PetscCall(PetscMalloc1(PETSC_MAX_PATH_LEN, tag_name));

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

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

    560:   /* increment the new value of n */
    561:   ++n;

    563: #ifdef MOAB_HAVE_MPI
    564:   /* Make sure that n is consistent across all processes */
    565:   PetscCallMPI(MPIU_Allreduce(&n, &global_n, 1, MPI_INT, MPI_MAX, pcomm->comm()));
    566: #else
    567:   global_n = n;
    568: #endif

    570:   /* Set the new name accordingly and return */
    571:   PetscCall(PetscSNPrintf(*tag_name, PETSC_MAX_PATH_LEN - 1, "%s_%" PetscInt_FMT, PVEC_PREFIX, global_n));
    572:   mberr = mbiface->tag_set_data(indexTag, &rootset, 1, (const void *)&global_n);
    573:   MBERRNM(mberr);
    574:   PetscFunctionReturn(PETSC_SUCCESS);
    575: }

    577: PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM dm, Vec *gvec)
    578: {
    579:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

    581:   PetscFunctionBegin;
    583:   PetscAssertPointer(gvec, 2);
    584:   PetscCall(DMCreateVector_Moab_Private(dm, NULL, dmmoab->vowned, PETSC_TRUE, PETSC_TRUE, gvec));
    585:   PetscFunctionReturn(PETSC_SUCCESS);
    586: }

    588: PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM dm, Vec *lvec)
    589: {
    590:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

    592:   PetscFunctionBegin;
    594:   PetscAssertPointer(lvec, 2);
    595:   PetscCall(DMCreateVector_Moab_Private(dm, NULL, dmmoab->vlocal, PETSC_FALSE, PETSC_TRUE, lvec));
    596:   PetscFunctionReturn(PETSC_SUCCESS);
    597: }

    599: static PetscErrorCode DMVecDuplicate_Moab(Vec x, Vec *y)
    600: {
    601:   DM             dm;
    602:   PetscContainer moabdata;
    603:   Vec_MOAB      *vmoab;

    605:   PetscFunctionBegin;
    607:   PetscAssertPointer(y, 2);

    609:   /* Get the Vec_MOAB struct for the original vector */
    610:   PetscCall(PetscObjectQuery((PetscObject)x, "MOABData", (PetscObject *)&moabdata));
    611:   PetscCall(PetscContainerGetPointer(moabdata, (void **)&vmoab));

    613:   PetscCall(VecGetDM(x, &dm));

    616:   PetscCall(DMCreateVector_Moab_Private(dm, NULL, vmoab->tag_range, vmoab->is_global_vec, PETSC_TRUE, y));
    617:   PetscCall(VecSetDM(*y, dm));
    618:   PetscFunctionReturn(PETSC_SUCCESS);
    619: }

    621: static PetscErrorCode DMVecUserDestroy_Moab(void *user)
    622: {
    623:   Vec_MOAB       *vmoab = (Vec_MOAB *)user;
    624:   moab::ErrorCode merr;

    626:   PetscFunctionBegin;
    627:   if (vmoab->new_tag && vmoab->tag) {
    628:     /* Tag was created via a call to VecDuplicate, delete the underlying tag in MOAB */
    629:     merr = vmoab->mbiface->tag_delete(vmoab->tag);
    630:     MBERRNM(merr);
    631:   }
    632:   delete vmoab->tag_range;
    633:   vmoab->tag     = NULL;
    634:   vmoab->mbiface = NULL;
    635: #ifdef MOAB_HAVE_MPI
    636:   vmoab->pcomm = NULL;
    637: #endif
    638:   PetscCall(PetscFree(vmoab));
    639:   PetscFunctionReturn(PETSC_SUCCESS);
    640: }

    642: PETSC_EXTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM dm, Vec g, InsertMode mode, Vec l)
    643: {
    644:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

    646:   PetscFunctionBegin;
    647:   PetscCall(VecScatterBegin(dmmoab->ltog_sendrecv, g, l, mode, SCATTER_REVERSE));
    648:   PetscFunctionReturn(PETSC_SUCCESS);
    649: }

    651: PETSC_EXTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM dm, Vec g, InsertMode mode, Vec l)
    652: {
    653:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

    655:   PetscFunctionBegin;
    656:   PetscCall(VecScatterEnd(dmmoab->ltog_sendrecv, g, l, mode, SCATTER_REVERSE));
    657:   PetscFunctionReturn(PETSC_SUCCESS);
    658: }

    660: PETSC_EXTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM dm, Vec l, InsertMode mode, Vec g)
    661: {
    662:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

    664:   PetscFunctionBegin;
    665:   PetscCall(VecScatterBegin(dmmoab->ltog_sendrecv, l, g, mode, SCATTER_FORWARD));
    666:   PetscFunctionReturn(PETSC_SUCCESS);
    667: }

    669: PETSC_EXTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM dm, Vec l, InsertMode mode, Vec g)
    670: {
    671:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

    673:   PetscFunctionBegin;
    674:   PetscCall(VecScatterEnd(dmmoab->ltog_sendrecv, l, g, mode, SCATTER_FORWARD));
    675:   PetscFunctionReturn(PETSC_SUCCESS);
    676: }