Actual source code: dmmbvec.cxx
petsc-3.10.5 2019-03-28
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: }