Actual source code: dmmbvec.cxx

petsc-3.11.4 2019-09-28
Report Typos and Errors
  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 on MPI_Comm

 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: .keywords: Vec, create

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

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

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


 52: /*@C
 53:   DMMoabGetVecTag - Get the MOAB tag associated with this Vec

 55:   Input Parameter:
 56: . vec - Vec being queried

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

 61:   Level: beginner

 63: .keywords: DMMoab, MOAB tag

 65: .seealso: DMMoabCreateVector(), DMMoabGetVecRange()
 66: @*/
 67: PetscErrorCode DMMoabGetVecTag(Vec vec, moab::Tag *tag)
 68: {
 69:   PetscContainer  moabdata;
 70:   Vec_MOAB        *vmoab;
 71:   PetscErrorCode  ierr;


 76:   /* Get the MOAB private data */
 77:   PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject*) &moabdata);
 78:   PetscContainerGetPointer(moabdata, (void**) &vmoab);

 80:   *tag = vmoab->tag;
 81:   return(0);
 82: }


 85: /*@C
 86:   DMMoabGetVecRange - Get the MOAB entities associated with this Vec

 88:   Input Parameter:
 89: . vec   - Vec being queried

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

 94:   Level: beginner

 96: .keywords: DMMoab, MOAB range

 98: .seealso: DMMoabCreateVector(), DMMoabGetVecTag()
 99: @*/
100: PetscErrorCode DMMoabGetVecRange(Vec vec, moab::Range *range)
101: {
102:   PetscContainer  moabdata;
103:   Vec_MOAB        *vmoab;
104:   PetscErrorCode  ierr;


109:   /* Get the MOAB private data handle */
110:   PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject*) &moabdata);
111:   PetscContainerGetPointer(moabdata, (void**) &vmoab);

113:   *range = *vmoab->tag_range;
114:   return(0);
115: }


118: /*@C
119:   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

121:   Collective on MPI_Comm

123:   Input Parameter:
124: . dm              - The DMMoab object being set
125: . vec             - The Vector whose underlying data is requested

127:   Output Parameter:
128: . array           - The local data array

130:   Level: intermediate

132: .keywords: MOAB, distributed array

134: .seealso: DMMoabVecRestoreArray(), DMMoabVecGetArrayRead(), DMMoabVecRestoreArrayRead()
135: @*/
136: PetscErrorCode  DMMoabVecGetArray(DM dm, Vec vec, void* array)
137: {
138:   DM_Moab        *dmmoab;
139:   moab::ErrorCode merr;
140:   PetscErrorCode  ierr;
141:   PetscInt        count, i, f;
142:   moab::Tag       vtag;
143:   PetscScalar     **varray;
144:   PetscScalar     *marray;
145:   PetscContainer  moabdata;
146:   Vec_MOAB        *vmoab, *xmoab;

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

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

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

161:   if (vmoab->is_native_vec) {

163:     /* get the local representation of the arrays from Vectors */
164:     DMGetLocalVector(dm, &vmoab->local);
165:     DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vmoab->local);
166:     DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vmoab->local);

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

172:     /* get the local representation of the arrays from Vectors */
173:     VecGhostGetLocalForm(vmoab->local, &xmoab->local);
174:     VecGhostUpdateBegin(vmoab->local, INSERT_VALUES, SCATTER_FORWARD);
175:     VecGhostUpdateEnd(vmoab->local, INSERT_VALUES, SCATTER_FORWARD);

177:     VecGetArray(xmoab->local, varray);
178:   }
179:   else {

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

184: #ifdef MOAB_HAVE_MPI
185:     /* exchange the data into ghost cells first */
186:     merr = dmmoab->pcomm->exchange_tags(vtag, *dmmoab->vlocal); MBERRNM(merr);
187: #endif

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

191:     /* Get the array data for local entities */
192:     merr = dmmoab->mbiface->tag_iterate(vtag, dmmoab->vlocal->begin(), dmmoab->vlocal->end(), count, reinterpret_cast<void*&>(marray), false); MBERRNM(merr);
193:     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);

195:     i = 0;
196:     for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
197:       for (f = 0; f < dmmoab->numFields; f++, i++)
198:         (*varray)[dmmoab->lidmap[(PetscInt)*iter - dmmoab->seqstart]*dmmoab->numFields + f] = marray[i];
199:       //(*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]]=marray[i];
200:     }
201:   }
202:   return(0);
203: }


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

209:   Collective on MPI_Comm

211:   Input Parameter:
212: + dm              - The DMMoab object being set
213: . vec             - The Vector whose underlying data is requested
214: . array           - The local data array

216:   Level: intermediate

218: .keywords: MOAB, distributed array

220: .seealso: DMMoabVecGetArray(), DMMoabVecGetArrayRead(), DMMoabVecRestoreArrayRead()
221: @*/
222: PetscErrorCode  DMMoabVecRestoreArray(DM dm, Vec vec, void* array)
223: {
224:   DM_Moab        *dmmoab;
225:   moab::ErrorCode merr;
226:   PetscErrorCode  ierr;
227:   moab::Tag       vtag;
228:   PetscInt        count, i, f;
229:   PetscScalar     **varray;
230:   PetscScalar     *marray;
231:   PetscContainer  moabdata;
232:   Vec_MOAB        *vmoab, *xmoab;

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

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

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

247:   if (vmoab->is_native_vec) {

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

253:     /* get the local representation of the arrays from Vectors */
254:     VecRestoreArray(xmoab->local, varray);
255:     VecGhostUpdateBegin(vmoab->local, ADD_VALUES, SCATTER_REVERSE);
256:     VecGhostUpdateEnd(vmoab->local, ADD_VALUES, SCATTER_REVERSE);
257:     VecGhostRestoreLocalForm(vmoab->local, &xmoab->local);

259:     /* restore local pieces */
260:     DMLocalToGlobalBegin(dm, vmoab->local, INSERT_VALUES, vec);
261:     DMLocalToGlobalEnd(dm, vmoab->local, INSERT_VALUES, vec);
262:     DMRestoreLocalVector(dm, &vmoab->local);
263:   }
264:   else {

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

269:     /* Get the array data for local entities */
270:     merr = dmmoab->mbiface->tag_iterate(vtag, dmmoab->vlocal->begin(), dmmoab->vlocal->end(), count, reinterpret_cast<void*&>(marray), false); MBERRNM(merr);
271:     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);

273:     i = 0;
274:     for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
275:       for (f = 0; f < dmmoab->numFields; f++, i++)
276:         marray[i] = (*varray)[dmmoab->lidmap[(PetscInt) * iter - dmmoab->seqstart] * dmmoab->numFields + f];
277:       //marray[i] = (*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]];
278:     }

280: #ifdef MOAB_HAVE_MPI
281:     /* reduce the tags correctly -> should probably let the user choose how to reduce in the future
282:       For all FEM residual based assembly calculations, MPI_SUM should serve well */
283:     merr = dmmoab->pcomm->reduce_tags(vtag, MPI_SUM, *dmmoab->vlocal); MBERRV(dmmoab->mbiface, merr);
284: #endif
285:     PetscFree(*varray);
286:   }
287:   return(0);
288: }

290: /*@C
291:   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

293:   Collective on MPI_Comm

295:   Input Parameter:
296: + dm              - The DMMoab object being set
297: . vec             - The Vector whose underlying data is requested

299:   Output Parameter:
300: . array           - The local data array

302:   Level: intermediate

304: .keywords: MOAB, distributed array

306: .seealso: DMMoabVecRestoreArrayRead(), DMMoabVecGetArray(), DMMoabVecRestoreArray()
307: @*/
308: PetscErrorCode  DMMoabVecGetArrayRead(DM dm, Vec vec, void* array)
309: {
310:   DM_Moab        *dmmoab;
311:   moab::ErrorCode merr;
312:   PetscErrorCode  ierr;
313:   PetscInt        count, i, f;
314:   moab::Tag       vtag;
315:   PetscScalar     **varray;
316:   PetscScalar     *marray;
317:   PetscContainer  moabdata;
318:   Vec_MOAB        *vmoab, *xmoab;

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

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

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

333:   if (vmoab->is_native_vec) {
334:     /* get the local representation of the arrays from Vectors */
335:     DMGetLocalVector(dm, &vmoab->local);
336:     DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vmoab->local);
337:     DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vmoab->local);

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

343:     /* get the local representation of the arrays from Vectors */
344:     VecGhostGetLocalForm(vmoab->local, &xmoab->local);
345:     VecGhostUpdateBegin(vmoab->local, INSERT_VALUES, SCATTER_FORWARD);
346:     VecGhostUpdateEnd(vmoab->local, INSERT_VALUES, SCATTER_FORWARD);
347:     VecGetArray(xmoab->local, varray);
348:   }
349:   else {
350:     /* Get the MOAB private data */
351:     DMMoabGetVecTag(vec, &vtag);

353: #ifdef MOAB_HAVE_MPI
354:     /* exchange the data into ghost cells first */
355:     merr = dmmoab->pcomm->exchange_tags(vtag, *dmmoab->vlocal); MBERRNM(merr);
356: #endif
357:     PetscMalloc1((dmmoab->nloc + dmmoab->nghost) * dmmoab->numFields, varray);

359:     /* Get the array data for local entities */
360:     merr = dmmoab->mbiface->tag_iterate(vtag, dmmoab->vlocal->begin(), dmmoab->vlocal->end(), count, reinterpret_cast<void*&>(marray), false); MBERRNM(merr);
361:     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);

363:     i = 0;
364:     for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
365:       for (f = 0; f < dmmoab->numFields; f++, i++)
366:         (*varray)[dmmoab->lidmap[(PetscInt)*iter - dmmoab->seqstart]*dmmoab->numFields + f] = marray[i];
367:       //(*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]]=marray[i];
368:     }
369:   }
370:   return(0);
371: }


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

377:   Collective on MPI_Comm

379:   Input Parameter:
380: + dm              - The DMMoab object being set
381: . vec             - The Vector whose underlying data is requested
382: . array           - The local data array

384:   Level: intermediate

386: .keywords: MOAB, distributed array

388: .seealso: DMMoabVecGetArrayRead(), DMMoabVecGetArray(), DMMoabVecRestoreArray()
389: @*/
390: PetscErrorCode  DMMoabVecRestoreArrayRead(DM dm, Vec vec, void* array)
391: {
392:   PetscErrorCode  ierr;
393:   PetscScalar     **varray;
394:   PetscContainer  moabdata;
395:   Vec_MOAB        *vmoab, *xmoab;


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

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

409:   if (vmoab->is_native_vec) {
410:     /* Get the Vec_MOAB struct for the original vector */
411:     PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject*) &moabdata);
412:     PetscContainerGetPointer(moabdata, (void**)&xmoab);

414:     /* restore the local representation of the arrays from Vectors */
415:     VecRestoreArray(xmoab->local, varray);
416:     VecGhostRestoreLocalForm(vmoab->local, &xmoab->local);

418:     /* restore local pieces */
419:     DMRestoreLocalVector(dm, &vmoab->local);
420:   }
421:   else {
422:     /* Nothing to do but just free the memory allocated before */
423:     PetscFree(*varray);

425:   }
426:   return(0);
427: }


430: PetscErrorCode DMCreateVector_Moab_Private(DM dm, moab::Tag tag, const moab::Range* userrange, PetscBool is_global_vec, PetscBool destroy_tag, Vec *vec)
431: {
432:   PetscErrorCode    ierr;
433:   moab::ErrorCode   merr;
434:   PetscBool         is_newtag;
435:   const moab::Range *range;
436:   PetscInt          count, lnative_vec, gnative_vec;
437:   std::string       ttname;
438:   PetscScalar       *data_ptr, *defaultvals;

440:   Vec_MOAB *vmoab;
441:   DM_Moab *dmmoab = (DM_Moab*)dm->data;
442: #ifdef MOAB_HAVE_MPI
443:   moab::ParallelComm *pcomm = dmmoab->pcomm;
444: #endif
445:   moab::Interface *mbiface = dmmoab->mbiface;

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

453: #ifndef USE_NATIVE_PETSCVEC
454:   /* If the tag data is in a single sequence, use PETSc native vector since tag_iterate isn't useful anymore */
455:   lnative_vec = (range->psize() - 1);
456: #else
457:   lnative_vec = 1; /* NOTE: Testing PETSc vector: will force to create native vector all the time */
458: //  lnative_vec=0; /* NOTE: Testing MOAB vector: will force to create MOAB tag_iterate based vector all the time */
459: #endif

461: #ifdef MOAB_HAVE_MPI
462:   MPIU_Allreduce(&lnative_vec, &gnative_vec, 1, MPI_INT, MPI_MAX, (((PetscObject)dm)->comm));
463: #else
464:   gnative_vec = lnative_vec;
465: #endif

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

471:   if (!vmoab->is_native_vec) {
472:     merr = moab::MB_SUCCESS;
473:     if (tag != 0) merr = mbiface->tag_get_name(tag, ttname);
474:     if (!ttname.length() || merr != moab::MB_SUCCESS) {
475:       /* get the new name for the anonymous MOABVec -> the tag_name will be destroyed along with Tag */
476:       char *tag_name = NULL;
477: #ifdef MOAB_HAVE_MPI
478:       DMVecCreateTagName_Moab_Private(mbiface,pcomm,&tag_name);
479: #else
480:       DMVecCreateTagName_Moab_Private(mbiface,&tag_name);
481: #endif
482:       is_newtag = PETSC_TRUE;

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

487:       /* Create the tag */
488:       merr = mbiface->tag_get_handle(tag_name, dmmoab->numFields, moab::MB_TYPE_DOUBLE, tag,
489:                                      moab::MB_TAG_DENSE | moab::MB_TAG_CREAT, defaultvals); MBERRNM(merr);
490:       PetscFree(tag_name);
491:       PetscFree(defaultvals);
492:     }
493:     else {
494:       /* Make sure the tag data is of type "double" */
495:       moab::DataType tag_type;
496:       merr = mbiface->tag_get_data_type(tag, tag_type); MBERRNM(merr);
497:       if (tag_type != moab::MB_TYPE_DOUBLE) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Tag data type must be MB_TYPE_DOUBLE");
498:       is_newtag = destroy_tag;
499:     }

501:     vmoab->tag = tag;
502:     vmoab->new_tag = is_newtag;
503:   }
504:   vmoab->mbiface = mbiface;
505: #ifdef MOAB_HAVE_MPI
506:   vmoab->pcomm = pcomm;
507: #endif
508:   vmoab->is_global_vec = is_global_vec;
509:   vmoab->tag_size = dmmoab->bs;

511:   if (vmoab->is_native_vec) {

513:     /* Create the PETSc Vector directly and attach our functions accordingly */
514:     if (!is_global_vec) {
515:       /* This is an MPI Vector with ghosted padding */
516:       VecCreateGhostBlock((((PetscObject)dm)->comm), dmmoab->bs, dmmoab->numFields * dmmoab->nloc,
517:                                  dmmoab->numFields * dmmoab->n, dmmoab->nghost, &dmmoab->gsindices[dmmoab->nloc], vec);
518:     }
519:     else {
520:       /* This is an MPI/SEQ Vector */
521:       VecCreate((((PetscObject)dm)->comm), vec);
522:       VecSetSizes(*vec, dmmoab->numFields * dmmoab->nloc, PETSC_DECIDE);
523:       VecSetBlockSize(*vec, dmmoab->bs);
524:       VecSetType(*vec, VECMPI);
525:     }
526:   }
527:   else {
528:     /* Call tag_iterate. This will cause MOAB to allocate memory for the
529:        tag data if it hasn't already happened */
530:     merr = mbiface->tag_iterate(tag, range->begin(), range->end(), count, (void*&)data_ptr); MBERRNM(merr);

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

536:     /* Create the PETSc Vector
537:       Query MOAB mesh to check if there are any ghosted entities
538:         -> if we do, create a ghosted vector to map correctly to the same layout
539:         -> else, create a non-ghosted parallel vector */
540:     if (!is_global_vec) {
541:       /* This is an MPI Vector with ghosted padding */
542:       VecCreateGhostBlockWithArray((((PetscObject)dm)->comm), dmmoab->bs, dmmoab->numFields * dmmoab->nloc,
543:                                           dmmoab->numFields * dmmoab->n, dmmoab->nghost, &dmmoab->gsindices[dmmoab->nloc], data_ptr, vec);
544:     }
545:     else {
546:       /* This is an MPI Vector without ghosted padding */
547:       VecCreateMPIWithArray((((PetscObject)dm)->comm), dmmoab->bs, dmmoab->numFields * range->size(),
548:                                    PETSC_DECIDE, data_ptr, vec);
549:     }
550:   }
551:   VecSetFromOptions(*vec);

553:   /* create a container and store the internal MOAB data for faster access based on Entities etc */
554:   PetscContainer moabdata;
555:   PetscContainerCreate(PETSC_COMM_WORLD, &moabdata);
556:   PetscContainerSetPointer(moabdata, vmoab);
557:   PetscContainerSetUserDestroy(moabdata, DMVecUserDestroy_Moab);
558:   PetscObjectCompose((PetscObject) * vec, "MOABData", (PetscObject)moabdata);
559:   (*vec)->ops->duplicate = DMVecDuplicate_Moab;
560:   PetscContainerDestroy(&moabdata);

562:   /* Vector created, manually set local to global mapping */
563:   if (dmmoab->ltog_map) {
564:     VecSetLocalToGlobalMapping(*vec, dmmoab->ltog_map);
565:   }

567:   /* set the DM reference to the vector */
568:   VecSetDM(*vec, dm);
569:   return(0);
570: }


573: /*  DMVecCreateTagName_Moab_Private
574:  *
575:  *  Creates a unique tag name that will be shared across processes. If
576:  *  pcomm is NULL, then this is a serial vector. A unique tag name
577:  *  will be returned in tag_name in either case.
578:  *
579:  *  The tag names have the format _PETSC_VEC_N where N is some integer.
580:  *
581:  *  NOTE: The tag_name is allocated in this routine; The user needs to free
582:  *        the character array.
583:  */
584: #ifdef MOAB_HAVE_MPI
585: PetscErrorCode DMVecCreateTagName_Moab_Private(moab::Interface *mbiface, moab::ParallelComm *pcomm, char** tag_name)
586: #else
587: PetscErrorCode DMVecCreateTagName_Moab_Private(moab::Interface *mbiface, char** tag_name)
588: #endif
589: {
590:   moab::ErrorCode mberr;
591:   PetscErrorCode  ierr;
592:   PetscInt        n, global_n;
593:   moab::Tag indexTag;

596:   const char*       PVEC_PREFIX      = "__PETSC_VEC_";
597:   PetscMalloc1(PETSC_MAX_PATH_LEN, tag_name);

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

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

609:   /* increment the new value of n */
610:   ++n;

612: #ifdef MOAB_HAVE_MPI
613:   /* Make sure that n is consistent across all processes */
614:   MPIU_Allreduce(&n, &global_n, 1, MPI_INT, MPI_MAX, pcomm->comm());
615: #else
616:   global_n = n;
617: #endif

619:   /* Set the new name accordingly and return */
620:   PetscSNPrintf(*tag_name, PETSC_MAX_PATH_LEN - 1, "%s_%D", PVEC_PREFIX, global_n);
621:   mberr = mbiface->tag_set_data(indexTag, &rootset, 1, (const void*)&global_n); MBERRNM(mberr);
622:   return(0);
623: }


626: PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM dm, Vec *gvec)
627: {
628:   PetscErrorCode  ierr;
629:   DM_Moab         *dmmoab = (DM_Moab*)dm->data;

634:   DMCreateVector_Moab_Private(dm,NULL,dmmoab->vowned,PETSC_TRUE,PETSC_TRUE,gvec);
635:   return(0);
636: }


639: PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM dm, Vec *lvec)
640: {
641:   PetscErrorCode  ierr;
642:   DM_Moab         *dmmoab = (DM_Moab*)dm->data;

647:   DMCreateVector_Moab_Private(dm,NULL,dmmoab->vlocal,PETSC_FALSE,PETSC_TRUE,lvec);
648:   return(0);
649: }


652: PetscErrorCode DMVecDuplicate_Moab(Vec x, Vec *y)
653: {
655:   DM             dm;
656:   PetscContainer  moabdata;
657:   Vec_MOAB        *vmoab;


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

667:   VecGetDM(x, &dm);

670:   DMCreateVector_Moab_Private(dm,NULL,vmoab->tag_range,vmoab->is_global_vec,PETSC_TRUE,y);
671:   VecSetDM(*y, dm);
672:   return(0);
673: }


676: PetscErrorCode DMVecUserDestroy_Moab(void *user)
677: {
678:   Vec_MOAB        *vmoab = (Vec_MOAB*)user;
679:   PetscErrorCode  ierr;
680:   moab::ErrorCode merr;

683:   if (vmoab->new_tag && vmoab->tag) {
684:     /* Tag was created via a call to VecDuplicate, delete the underlying tag in MOAB */
685:     merr = vmoab->mbiface->tag_delete(vmoab->tag); MBERRNM(merr);
686:   }
687:   delete vmoab->tag_range;
688:   vmoab->tag = NULL;
689:   vmoab->mbiface = NULL;
690: #ifdef MOAB_HAVE_MPI
691:   vmoab->pcomm = NULL;
692: #endif
693:   PetscFree(vmoab);
694:   return(0);
695: }


698: PETSC_EXTERN PetscErrorCode  DMGlobalToLocalBegin_Moab(DM dm, Vec g, InsertMode mode, Vec l)
699: {
700:   PetscErrorCode    ierr;
701:   DM_Moab         *dmmoab = (DM_Moab*)dm->data;

704:   VecScatterBegin(dmmoab->ltog_sendrecv, g, l, mode, SCATTER_REVERSE);
705:   return(0);
706: }


709: PETSC_EXTERN PetscErrorCode  DMGlobalToLocalEnd_Moab(DM dm, Vec g, InsertMode mode, Vec l)
710: {
711:   PetscErrorCode    ierr;
712:   DM_Moab         *dmmoab = (DM_Moab*)dm->data;

715:   VecScatterEnd(dmmoab->ltog_sendrecv, g, l, mode, SCATTER_REVERSE);
716:   return(0);
717: }


720: PETSC_EXTERN PetscErrorCode  DMLocalToGlobalBegin_Moab(DM dm, Vec l, InsertMode mode, Vec g)
721: {
722:   PetscErrorCode    ierr;
723:   DM_Moab         *dmmoab = (DM_Moab*)dm->data;

726:   VecScatterBegin(dmmoab->ltog_sendrecv, l, g, mode, SCATTER_FORWARD);
727:   return(0);
728: }


731: PETSC_EXTERN PetscErrorCode  DMLocalToGlobalEnd_Moab(DM dm, Vec l, InsertMode mode, Vec g)
732: {
733:   PetscErrorCode    ierr;
734:   DM_Moab         *dmmoab = (DM_Moab*)dm->data;

737:   VecScatterEnd(dmmoab->ltog_sendrecv, l, g, mode, SCATTER_FORWARD);
738:   return(0);
739: }