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:   PetscCall(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:   PetscContainer moabdata;
515:   PetscCall(PetscContainerCreate(PETSC_COMM_WORLD, &moabdata));
516:   PetscCall(PetscContainerSetPointer(moabdata, vmoab));
517:   PetscCall(PetscContainerSetUserDestroy(moabdata, DMVecUserDestroy_Moab));
518:   PetscCall(PetscObjectCompose((PetscObject)*vec, "MOABData", (PetscObject)moabdata));
519:   (*vec)->ops->duplicate = DMVecDuplicate_Moab;
520:   PetscCall(PetscContainerDestroy(&moabdata));

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

525:   /* set the DM reference to the vector */
526:   PetscCall(VecSetDM(*vec, dm));
527:   PetscFunctionReturn(PETSC_SUCCESS);
528: }

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

551:   PetscFunctionBegin;
552:   const char *PVEC_PREFIX = "__PETSC_VEC_";
553:   PetscCall(PetscMalloc1(PETSC_MAX_PATH_LEN, tag_name));

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

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

565:   /* increment the new value of n */
566:   ++n;

568: #ifdef MOAB_HAVE_MPI
569:   /* Make sure that n is consistent across all processes */
570:   PetscCall(MPIU_Allreduce(&n, &global_n, 1, MPI_INT, MPI_MAX, pcomm->comm()));
571: #else
572:   global_n = n;
573: #endif

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

582: PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM dm, Vec *gvec)
583: {
584:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

586:   PetscFunctionBegin;
588:   PetscAssertPointer(gvec, 2);
589:   PetscCall(DMCreateVector_Moab_Private(dm, NULL, dmmoab->vowned, PETSC_TRUE, PETSC_TRUE, gvec));
590:   PetscFunctionReturn(PETSC_SUCCESS);
591: }

593: PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM dm, Vec *lvec)
594: {
595:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

597:   PetscFunctionBegin;
599:   PetscAssertPointer(lvec, 2);
600:   PetscCall(DMCreateVector_Moab_Private(dm, NULL, dmmoab->vlocal, PETSC_FALSE, PETSC_TRUE, lvec));
601:   PetscFunctionReturn(PETSC_SUCCESS);
602: }

604: static PetscErrorCode DMVecDuplicate_Moab(Vec x, Vec *y)
605: {
606:   DM             dm;
607:   PetscContainer moabdata;
608:   Vec_MOAB      *vmoab;

610:   PetscFunctionBegin;
612:   PetscAssertPointer(y, 2);

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

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

621:   PetscCall(DMCreateVector_Moab_Private(dm, NULL, vmoab->tag_range, vmoab->is_global_vec, PETSC_TRUE, y));
622:   PetscCall(VecSetDM(*y, dm));
623:   PetscFunctionReturn(PETSC_SUCCESS);
624: }

626: static PetscErrorCode DMVecUserDestroy_Moab(void *user)
627: {
628:   Vec_MOAB       *vmoab = (Vec_MOAB *)user;
629:   moab::ErrorCode merr;

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

647: PETSC_EXTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM dm, Vec g, InsertMode mode, Vec l)
648: {
649:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

651:   PetscFunctionBegin;
652:   PetscCall(VecScatterBegin(dmmoab->ltog_sendrecv, g, l, mode, SCATTER_REVERSE));
653:   PetscFunctionReturn(PETSC_SUCCESS);
654: }

656: PETSC_EXTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM dm, Vec g, InsertMode mode, Vec l)
657: {
658:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

660:   PetscFunctionBegin;
661:   PetscCall(VecScatterEnd(dmmoab->ltog_sendrecv, g, l, mode, SCATTER_REVERSE));
662:   PetscFunctionReturn(PETSC_SUCCESS);
663: }

665: PETSC_EXTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM dm, Vec l, InsertMode mode, Vec g)
666: {
667:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

669:   PetscFunctionBegin;
670:   PetscCall(VecScatterBegin(dmmoab->ltog_sendrecv, l, g, mode, SCATTER_FORWARD));
671:   PetscFunctionReturn(PETSC_SUCCESS);
672: }

674: PETSC_EXTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM dm, Vec l, InsertMode mode, Vec g)
675: {
676:   DM_Moab *dmmoab = (DM_Moab *)dm->data;

678:   PetscFunctionBegin;
679:   PetscCall(VecScatterEnd(dmmoab->ltog_sendrecv, l, g, mode, SCATTER_FORWARD));
680:   PetscFunctionReturn(PETSC_SUCCESS);
681: }