Actual source code: dmmbvec.cxx

petsc-3.14.6 2021-03-30
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

 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: }