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