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, moab::Tag, const moab::Range *, PetscBool, PetscBool, Vec *);
11: static PetscErrorCode DMVecCtxDestroy_Moab(PetscCtxRt);
12: static PetscErrorCode DMVecDuplicate_Moab(Vec, Vec *);
13: #ifdef MOAB_HAVE_MPI
14: static PetscErrorCode DMVecCreateTagName_Moab_Private(moab::Interface *, moab::ParallelComm *, char **);
15: #else
16: static PetscErrorCode DMVecCreateTagName_Moab_Private(moab::Interface *, char **);
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, &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, &vmoab));
100: *range = *vmoab->tag_range;
101: PetscFunctionReturn(PETSC_SUCCESS);
102: }
104: /*@C
105: 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
107: Collective
109: Input Parameters:
110: + dm - The DMMoab object being set
111: - vec - The Vector whose underlying data is requested
113: Output Parameter:
114: . array - The local data array
116: Level: intermediate
118: .seealso: `DMMoabVecRestoreArray()`, `DMMoabVecGetArrayRead()`, `DMMoabVecRestoreArrayRead()`
119: @*/
120: PetscErrorCode DMMoabVecGetArray(DM dm, Vec vec, void *array)
121: {
122: DM_Moab *dmmoab;
123: PetscInt count, i, f;
124: moab::Tag vtag;
125: PetscScalar **varray;
126: PetscScalar *marray;
127: PetscContainer moabdata;
128: Vec_MOAB *vmoab, *xmoab;
130: PetscFunctionBegin;
133: PetscAssertPointer(array, 3);
134: dmmoab = (DM_Moab *)dm->data;
136: /* Get the Vec_MOAB struct for the original vector */
137: PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
138: PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
140: /* Get the real scalar array handle */
141: varray = reinterpret_cast<PetscScalar **>(array);
143: if (vmoab->is_native_vec) {
144: /* get the local representation of the arrays from Vectors */
145: PetscCall(DMGetLocalVector(dm, &vmoab->local));
146: PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vmoab->local));
147: PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vmoab->local));
149: /* Get the Vec_MOAB struct for the original vector */
150: PetscCall(PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject *)&moabdata));
151: PetscCall(PetscContainerGetPointer(moabdata, &xmoab));
153: /* get the local representation of the arrays from Vectors */
154: PetscCall(VecGhostGetLocalForm(vmoab->local, &xmoab->local));
155: PetscCall(VecGhostUpdateBegin(vmoab->local, INSERT_VALUES, SCATTER_FORWARD));
156: PetscCall(VecGhostUpdateEnd(vmoab->local, INSERT_VALUES, SCATTER_FORWARD));
158: PetscCall(VecGetArray(xmoab->local, varray));
159: } else {
160: /* Get the MOAB private data */
161: PetscCall(DMMoabGetVecTag(vec, &vtag));
163: #ifdef MOAB_HAVE_MPI
164: /* exchange the data into ghost cells first */
165: PetscCallMOAB(dmmoab->pcomm->exchange_tags(vtag, *dmmoab->vlocal));
166: #endif
168: PetscCall(PetscMalloc1((dmmoab->nloc + dmmoab->nghost) * dmmoab->numFields, varray));
170: /* Get the array data for local entities */
171: PetscCallMOAB(dmmoab->mbiface->tag_iterate(vtag, dmmoab->vlocal->begin(), dmmoab->vlocal->end(), count, reinterpret_cast<void *&>(marray), false));
172: 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);
174: i = 0;
175: for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
176: for (f = 0; f < dmmoab->numFields; f++, i++) (*varray)[dmmoab->lidmap[(PetscInt)*iter - dmmoab->seqstart] * dmmoab->numFields + f] = marray[i];
177: //(*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]]=marray[i];
178: }
179: }
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: /*@C
184: DMMoabVecRestoreArray - Restores the writable direct access array obtained via DMMoabVecGetArray
186: Collective
188: Input Parameters:
189: + dm - The DMMoab object being set
190: . vec - The Vector whose underlying data is requested
191: - array - The local data array
193: Level: intermediate
195: .seealso: `DMMoabVecGetArray()`, `DMMoabVecGetArrayRead()`, `DMMoabVecRestoreArrayRead()`
196: @*/
197: PetscErrorCode DMMoabVecRestoreArray(DM dm, Vec vec, void *array)
198: {
199: DM_Moab *dmmoab;
200: moab::Tag vtag;
201: PetscInt count, i, f;
202: PetscScalar **varray;
203: PetscScalar *marray;
204: PetscContainer moabdata;
205: Vec_MOAB *vmoab, *xmoab;
207: PetscFunctionBegin;
210: PetscAssertPointer(array, 3);
211: dmmoab = (DM_Moab *)dm->data;
213: /* Get the Vec_MOAB struct for the original vector */
214: PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
215: PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
217: /* Get the real scalar array handle */
218: varray = reinterpret_cast<PetscScalar **>(array);
220: if (vmoab->is_native_vec) {
221: /* Get the Vec_MOAB struct for the original vector */
222: PetscCall(PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject *)&moabdata));
223: PetscCall(PetscContainerGetPointer(moabdata, &xmoab));
225: /* get the local representation of the arrays from Vectors */
226: PetscCall(VecRestoreArray(xmoab->local, varray));
227: PetscCall(VecGhostUpdateBegin(vmoab->local, ADD_VALUES, SCATTER_REVERSE));
228: PetscCall(VecGhostUpdateEnd(vmoab->local, ADD_VALUES, SCATTER_REVERSE));
229: PetscCall(VecGhostRestoreLocalForm(vmoab->local, &xmoab->local));
231: /* restore local pieces */
232: PetscCall(DMLocalToGlobalBegin(dm, vmoab->local, INSERT_VALUES, vec));
233: PetscCall(DMLocalToGlobalEnd(dm, vmoab->local, INSERT_VALUES, vec));
234: PetscCall(DMRestoreLocalVector(dm, &vmoab->local));
235: } else {
236: /* Get the MOAB private data */
237: PetscCall(DMMoabGetVecTag(vec, &vtag));
239: /* Get the array data for local entities */
240: PetscCallMOAB(dmmoab->mbiface->tag_iterate(vtag, dmmoab->vlocal->begin(), dmmoab->vlocal->end(), count, reinterpret_cast<void *&>(marray), false));
241: 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);
243: i = 0;
244: for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
245: for (f = 0; f < dmmoab->numFields; f++, i++) marray[i] = (*varray)[dmmoab->lidmap[(PetscInt)*iter - dmmoab->seqstart] * dmmoab->numFields + f];
246: //marray[i] = (*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]];
247: }
249: #ifdef MOAB_HAVE_MPI
250: /* reduce the tags correctly -> should probably let the user choose how to reduce in the future
251: For all FEM residual based assembly calculations, MPI_SUM should serve well */
252: PetscCallMOAB(dmmoab->pcomm->reduce_tags(vtag, MPI_SUM, *dmmoab->vlocal));
253: #endif
254: PetscCall(PetscFree(*varray));
255: }
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: /*@C
260: 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
262: Collective
264: Input Parameters:
265: + dm - The DMMoab object being set
266: - vec - The Vector whose underlying data is requested
268: Output Parameter:
269: . array - The local data array
271: Level: intermediate
273: .seealso: `DMMoabVecRestoreArrayRead()`, `DMMoabVecGetArray()`, `DMMoabVecRestoreArray()`
274: @*/
275: PetscErrorCode DMMoabVecGetArrayRead(DM dm, Vec vec, void *array)
276: {
277: DM_Moab *dmmoab;
278: PetscInt count, i, f;
279: moab::Tag vtag;
280: PetscScalar **varray;
281: PetscScalar *marray;
282: PetscContainer moabdata;
283: Vec_MOAB *vmoab, *xmoab;
285: PetscFunctionBegin;
288: PetscAssertPointer(array, 3);
289: dmmoab = (DM_Moab *)dm->data;
291: /* Get the Vec_MOAB struct for the original vector */
292: PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
293: PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
295: /* Get the real scalar array handle */
296: varray = reinterpret_cast<PetscScalar **>(array);
298: if (vmoab->is_native_vec) {
299: /* get the local representation of the arrays from Vectors */
300: PetscCall(DMGetLocalVector(dm, &vmoab->local));
301: PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vmoab->local));
302: PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vmoab->local));
304: /* Get the Vec_MOAB struct for the original vector */
305: PetscCall(PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject *)&moabdata));
306: PetscCall(PetscContainerGetPointer(moabdata, &xmoab));
308: /* get the local representation of the arrays from Vectors */
309: PetscCall(VecGhostGetLocalForm(vmoab->local, &xmoab->local));
310: PetscCall(VecGhostUpdateBegin(vmoab->local, INSERT_VALUES, SCATTER_FORWARD));
311: PetscCall(VecGhostUpdateEnd(vmoab->local, INSERT_VALUES, SCATTER_FORWARD));
312: PetscCall(VecGetArray(xmoab->local, varray));
313: } else {
314: /* Get the MOAB private data */
315: PetscCall(DMMoabGetVecTag(vec, &vtag));
317: #ifdef MOAB_HAVE_MPI
318: /* exchange the data into ghost cells first */
319: PetscCallMOAB(dmmoab->pcomm->exchange_tags(vtag, *dmmoab->vlocal));
320: #endif
321: PetscCall(PetscMalloc1((dmmoab->nloc + dmmoab->nghost) * dmmoab->numFields, varray));
323: /* Get the array data for local entities */
324: PetscCallMOAB(dmmoab->mbiface->tag_iterate(vtag, dmmoab->vlocal->begin(), dmmoab->vlocal->end(), count, reinterpret_cast<void *&>(marray), false));
325: 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);
327: i = 0;
328: for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
329: for (f = 0; f < dmmoab->numFields; f++, i++) (*varray)[dmmoab->lidmap[(PetscInt)*iter - dmmoab->seqstart] * dmmoab->numFields + f] = marray[i];
330: //(*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]]=marray[i];
331: }
332: }
333: PetscFunctionReturn(PETSC_SUCCESS);
334: }
336: /*@C
337: DMMoabVecRestoreArrayRead - Restores the read-only direct access array obtained via DMMoabVecGetArray
339: Collective
341: Input Parameters:
342: + dm - The DMMoab object being set
343: . vec - The Vector whose underlying data is requested
344: - array - The local data array
346: Level: intermediate
348: .seealso: `DMMoabVecGetArrayRead()`, `DMMoabVecGetArray()`, `DMMoabVecRestoreArray()`
349: @*/
350: PetscErrorCode DMMoabVecRestoreArrayRead(DM dm, Vec vec, void *array)
351: {
352: PetscScalar **varray;
353: PetscContainer moabdata;
354: Vec_MOAB *vmoab, *xmoab;
356: PetscFunctionBegin;
359: PetscAssertPointer(array, 3);
361: /* Get the Vec_MOAB struct for the original vector */
362: PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
363: PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
365: /* Get the real scalar array handle */
366: varray = reinterpret_cast<PetscScalar **>(array);
368: if (vmoab->is_native_vec) {
369: /* Get the Vec_MOAB struct for the original vector */
370: PetscCall(PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject *)&moabdata));
371: PetscCall(PetscContainerGetPointer(moabdata, &xmoab));
373: /* restore the local representation of the arrays from Vectors */
374: PetscCall(VecRestoreArray(xmoab->local, varray));
375: PetscCall(VecGhostRestoreLocalForm(vmoab->local, &xmoab->local));
377: /* restore local pieces */
378: PetscCall(DMRestoreLocalVector(dm, &vmoab->local));
379: } else {
380: /* Nothing to do but just free the memory allocated before */
381: PetscCall(PetscFree(*varray));
382: }
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: static PetscErrorCode DMCreateVector_Moab_Private(DM dm, moab::Tag tag, const moab::Range *userrange, PetscBool is_global_vec, PetscBool destroy_tag, Vec *vec)
387: {
388: moab::ErrorCode merr;
389: PetscBool is_newtag;
390: const moab::Range *range;
391: PetscInt count, lnative_vec, gnative_vec;
392: std::string ttname;
393: PetscScalar *data_ptr, *defaultvals;
395: Vec_MOAB *vmoab;
396: DM_Moab *dmmoab = (DM_Moab *)dm->data;
397: #ifdef MOAB_HAVE_MPI
398: moab::ParallelComm *pcomm = dmmoab->pcomm;
399: #endif
400: moab::Interface *mbiface = dmmoab->mbiface;
402: PetscFunctionBegin;
403: PetscCheck(sizeof(PetscReal) == sizeof(PetscScalar), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "MOAB tags only support Real types (Complex-type unsupported)");
404: if (!userrange) range = dmmoab->vowned;
405: else range = userrange;
406: PetscCheck(range, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Input range cannot be empty or call DMSetUp first.");
408: #if !defined(USE_NATIVE_PETSCVEC)
409: /* If the tag data is in a single sequence, use PETSc native vector since tag_iterate isn't useful anymore */
410: lnative_vec = (range->psize() - 1);
411: #else
412: lnative_vec = 1; /* NOTE: Testing PETSc vector: will force to create native vector all the time */
413: // lnative_vec=0; /* NOTE: Testing MOAB vector: will force to create MOAB tag_iterate based vector all the time */
414: #endif
416: #ifdef MOAB_HAVE_MPI
417: PetscCallMPI(MPIU_Allreduce(&lnative_vec, &gnative_vec, 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm));
418: #else
419: gnative_vec = lnative_vec;
420: #endif
422: /* Create the MOAB internal data object */
423: PetscCall(PetscNew(&vmoab));
424: vmoab->is_native_vec = (gnative_vec > 0 ? PETSC_TRUE : PETSC_FALSE);
426: if (!vmoab->is_native_vec) {
427: merr = moab::MB_SUCCESS;
428: if (tag != 0) merr = mbiface->tag_get_name(tag, ttname);
429: if (!ttname.length() || merr != moab::MB_SUCCESS) {
430: /* get the new name for the anonymous MOABVec -> the tag_name will be destroyed along with Tag */
431: char *tag_name = NULL;
432: #ifdef MOAB_HAVE_MPI
433: PetscCall(DMVecCreateTagName_Moab_Private(mbiface, pcomm, &tag_name));
434: #else
435: PetscCall(DMVecCreateTagName_Moab_Private(mbiface, &tag_name));
436: #endif
437: is_newtag = PETSC_TRUE;
439: /* Create the default value for the tag (all zeros) */
440: PetscCall(PetscCalloc1(dmmoab->numFields, &defaultvals));
442: /* Create the tag */
443: PetscCallMOAB(mbiface->tag_get_handle(tag_name, dmmoab->numFields, moab::MB_TYPE_DOUBLE, tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT, defaultvals));
444: PetscCall(PetscFree(tag_name));
445: PetscCall(PetscFree(defaultvals));
446: } else {
447: /* Make sure the tag data is of type "double" */
448: moab::DataType tag_type;
449: PetscCallMOAB(mbiface->tag_get_data_type(tag, tag_type));
450: PetscCheck(tag_type == moab::MB_TYPE_DOUBLE, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Tag data type must be MB_TYPE_DOUBLE");
451: is_newtag = destroy_tag;
452: }
454: vmoab->tag = tag;
455: vmoab->new_tag = is_newtag;
456: }
457: vmoab->mbiface = mbiface;
458: #ifdef MOAB_HAVE_MPI
459: vmoab->pcomm = pcomm;
460: #endif
461: vmoab->is_global_vec = is_global_vec;
462: vmoab->tag_size = dmmoab->bs;
464: if (vmoab->is_native_vec) {
465: /* Create the PETSc Vector directly and attach our functions accordingly */
466: if (!is_global_vec) {
467: /* This is an MPI Vector with ghosted padding */
468: PetscCall(VecCreateGhostBlock(((PetscObject)dm)->comm, dmmoab->bs, dmmoab->numFields * dmmoab->nloc, dmmoab->numFields * dmmoab->n, dmmoab->nghost, &dmmoab->gsindices[dmmoab->nloc], vec));
469: } else {
470: /* This is an MPI/SEQ Vector */
471: PetscCall(VecCreate(((PetscObject)dm)->comm, vec));
472: PetscCall(VecSetSizes(*vec, dmmoab->numFields * dmmoab->nloc, PETSC_DECIDE));
473: PetscCall(VecSetBlockSize(*vec, dmmoab->bs));
474: PetscCall(VecSetType(*vec, VECMPI));
475: }
476: } else {
477: /* Call tag_iterate. This will cause MOAB to allocate memory for the
478: tag data if it hasn't already happened */
479: PetscCallMOAB(mbiface->tag_iterate(tag, range->begin(), range->end(), count, (void *&)data_ptr));
481: /* set the reference for vector range */
482: vmoab->tag_range = new moab::Range(*range);
483: PetscCallMOAB(mbiface->tag_get_length(tag, dmmoab->numFields));
485: /* Create the PETSc Vector
486: Query MOAB mesh to check if there are any ghosted entities
487: -> if we do, create a ghosted vector to map correctly to the same layout
488: -> else, create a non-ghosted parallel vector */
489: if (!is_global_vec) {
490: /* This is an MPI Vector with ghosted padding */
491: PetscCall(VecCreateGhostBlockWithArray(((PetscObject)dm)->comm, dmmoab->bs, dmmoab->numFields * dmmoab->nloc, dmmoab->numFields * dmmoab->n, dmmoab->nghost, &dmmoab->gsindices[dmmoab->nloc], data_ptr, vec));
492: } else {
493: /* This is an MPI Vector without ghosted padding */
494: PetscCall(VecCreateMPIWithArray(((PetscObject)dm)->comm, dmmoab->bs, dmmoab->numFields * range->size(), PETSC_DECIDE, data_ptr, vec));
495: }
496: }
497: PetscCall(VecSetFromOptions(*vec));
499: /* create a container and store the internal MOAB data for faster access based on Entities etc */
500: PetscCall(PetscObjectContainerCompose((PetscObject)*vec, "MOABData", vmoab, DMVecCtxDestroy_Moab));
501: (*vec)->ops->duplicate = DMVecDuplicate_Moab;
503: /* Vector created, manually set local to global mapping */
504: if (dmmoab->ltog_map) PetscCall(VecSetLocalToGlobalMapping(*vec, dmmoab->ltog_map));
506: /* set the DM reference to the vector */
507: PetscCall(VecSetDM(*vec, dm));
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: /* DMVecCreateTagName_Moab_Private
512: *
513: * Creates a unique tag name that will be shared across processes. If
514: * pcomm is NULL, then this is a serial vector. A unique tag name
515: * will be returned in tag_name in either case.
516: *
517: * The tag names have the format _PETSC_VEC_N where N is some integer.
518: *
519: * NOTE: The tag_name is allocated in this routine; The user needs to free
520: * the character array.
521: */
522: #ifdef MOAB_HAVE_MPI
523: static PetscErrorCode DMVecCreateTagName_Moab_Private(moab::Interface *mbiface, moab::ParallelComm *pcomm, char **tag_name)
524: #else
525: static PetscErrorCode DMVecCreateTagName_Moab_Private(moab::Interface *mbiface, char **tag_name)
526: #endif
527: {
528: moab::ErrorCode mberr;
529: PetscInt n, global_n;
530: moab::Tag indexTag;
532: PetscFunctionBegin;
533: const char *PVEC_PREFIX = "__PETSC_VEC_";
534: PetscCall(PetscMalloc1(PETSC_MAX_PATH_LEN, tag_name));
536: moab::EntityHandle rootset = mbiface->get_root_set();
538: /* Check to see if there are any PETSc vectors defined */
539: /* Create a tag in MOAB mesh to index and keep track of number of PETSc vec tags */
540: PetscCallMOAB(mbiface->tag_get_handle("__PETSC_VECS__", 1, moab::MB_TYPE_INTEGER, indexTag, moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT, 0));
541: mberr = mbiface->tag_get_data(indexTag, &rootset, 1, &n);
542: if (mberr == moab::MB_TAG_NOT_FOUND) n = 0; /* this is the first temporary vector */
543: else PetscCheck(mberr == moab::MB_SUCCESS, PETSC_COMM_SELF, PETSC_ERR_LIB, "MOAB error %d in tag_get_data", mberr);
545: /* increment the new value of n */
546: ++n;
548: #ifdef MOAB_HAVE_MPI
549: /* Make sure that n is consistent across all processes */
550: PetscCallMPI(MPIU_Allreduce(&n, &global_n, 1, MPI_INT, MPI_MAX, pcomm->comm()));
551: #else
552: global_n = n;
553: #endif
555: /* Set the new name accordingly and return */
556: PetscCall(PetscSNPrintf(*tag_name, PETSC_MAX_PATH_LEN - 1, "%s_%" PetscInt_FMT, PVEC_PREFIX, global_n));
557: PetscCallMOAB(mbiface->tag_set_data(indexTag, &rootset, 1, (const void *)&global_n));
558: PetscFunctionReturn(PETSC_SUCCESS);
559: }
561: PETSC_INTERN PetscErrorCode DMCreateGlobalVector_Moab(DM dm, Vec *gvec)
562: {
563: DM_Moab *dmmoab = (DM_Moab *)dm->data;
565: PetscFunctionBegin;
567: PetscAssertPointer(gvec, 2);
568: PetscCall(DMCreateVector_Moab_Private(dm, NULL, dmmoab->vowned, PETSC_TRUE, PETSC_TRUE, gvec));
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: PETSC_INTERN PetscErrorCode DMCreateLocalVector_Moab(DM dm, Vec *lvec)
573: {
574: DM_Moab *dmmoab = (DM_Moab *)dm->data;
576: PetscFunctionBegin;
578: PetscAssertPointer(lvec, 2);
579: PetscCall(DMCreateVector_Moab_Private(dm, NULL, dmmoab->vlocal, PETSC_FALSE, PETSC_TRUE, lvec));
580: PetscFunctionReturn(PETSC_SUCCESS);
581: }
583: static PetscErrorCode DMVecDuplicate_Moab(Vec x, Vec *y)
584: {
585: DM dm;
586: PetscContainer moabdata;
587: Vec_MOAB *vmoab;
589: PetscFunctionBegin;
591: PetscAssertPointer(y, 2);
593: /* Get the Vec_MOAB struct for the original vector */
594: PetscCall(PetscObjectQuery((PetscObject)x, "MOABData", (PetscObject *)&moabdata));
595: PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
597: PetscCall(VecGetDM(x, &dm));
600: PetscCall(DMCreateVector_Moab_Private(dm, NULL, vmoab->tag_range, vmoab->is_global_vec, PETSC_TRUE, y));
601: PetscCall(VecSetDM(*y, dm));
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: static PetscErrorCode DMVecCtxDestroy_Moab(PetscCtxRt ctx)
606: {
607: Vec_MOAB *vmoab = *(Vec_MOAB **)ctx;
609: PetscFunctionBegin;
610: /* Tag was created via a call to VecDuplicate, delete the underlying tag in MOAB */
611: if (vmoab->new_tag && vmoab->tag) PetscCallMOAB(vmoab->mbiface->tag_delete(vmoab->tag));
612: delete vmoab->tag_range;
613: vmoab->tag = NULL;
614: vmoab->mbiface = NULL;
615: #ifdef MOAB_HAVE_MPI
616: vmoab->pcomm = NULL;
617: #endif
618: PetscCall(PetscFree(vmoab));
619: PetscFunctionReturn(PETSC_SUCCESS);
620: }
622: PETSC_INTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM dm, Vec g, InsertMode mode, Vec l)
623: {
624: DM_Moab *dmmoab = (DM_Moab *)dm->data;
626: PetscFunctionBegin;
627: PetscCall(VecScatterBegin(dmmoab->ltog_sendrecv, g, l, mode, SCATTER_REVERSE));
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: PETSC_INTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM dm, Vec g, InsertMode mode, Vec l)
632: {
633: DM_Moab *dmmoab = (DM_Moab *)dm->data;
635: PetscFunctionBegin;
636: PetscCall(VecScatterEnd(dmmoab->ltog_sendrecv, g, l, mode, SCATTER_REVERSE));
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }
640: PETSC_INTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM dm, Vec l, InsertMode mode, Vec g)
641: {
642: DM_Moab *dmmoab = (DM_Moab *)dm->data;
644: PetscFunctionBegin;
645: PetscCall(VecScatterBegin(dmmoab->ltog_sendrecv, l, g, mode, SCATTER_FORWARD));
646: PetscFunctionReturn(PETSC_SUCCESS);
647: }
649: PETSC_INTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM dm, Vec l, InsertMode mode, Vec g)
650: {
651: DM_Moab *dmmoab = (DM_Moab *)dm->data;
653: PetscFunctionBegin;
654: PetscCall(VecScatterEnd(dmmoab->ltog_sendrecv, l, g, mode, SCATTER_FORWARD));
655: PetscFunctionReturn(PETSC_SUCCESS);
656: }