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: moab::ErrorCode merr;
124: PetscInt count, i, f;
125: moab::Tag vtag;
126: PetscScalar **varray;
127: PetscScalar *marray;
128: PetscContainer moabdata;
129: Vec_MOAB *vmoab, *xmoab;
131: PetscFunctionBegin;
134: PetscAssertPointer(array, 3);
135: dmmoab = (DM_Moab *)dm->data;
137: /* Get the Vec_MOAB struct for the original vector */
138: PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
139: PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
141: /* Get the real scalar array handle */
142: varray = reinterpret_cast<PetscScalar **>(array);
144: if (vmoab->is_native_vec) {
145: /* get the local representation of the arrays from Vectors */
146: PetscCall(DMGetLocalVector(dm, &vmoab->local));
147: PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vmoab->local));
148: PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vmoab->local));
150: /* Get the Vec_MOAB struct for the original vector */
151: PetscCall(PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject *)&moabdata));
152: PetscCall(PetscContainerGetPointer(moabdata, &xmoab));
154: /* get the local representation of the arrays from Vectors */
155: PetscCall(VecGhostGetLocalForm(vmoab->local, &xmoab->local));
156: PetscCall(VecGhostUpdateBegin(vmoab->local, INSERT_VALUES, SCATTER_FORWARD));
157: PetscCall(VecGhostUpdateEnd(vmoab->local, INSERT_VALUES, SCATTER_FORWARD));
159: PetscCall(VecGetArray(xmoab->local, varray));
160: } else {
161: /* Get the MOAB private data */
162: PetscCall(DMMoabGetVecTag(vec, &vtag));
164: #ifdef MOAB_HAVE_MPI
165: /* exchange the data into ghost cells first */
166: merr = dmmoab->pcomm->exchange_tags(vtag, *dmmoab->vlocal);
167: MBERRNM(merr);
168: #endif
170: PetscCall(PetscMalloc1((dmmoab->nloc + dmmoab->nghost) * dmmoab->numFields, varray));
172: /* Get the array data for local entities */
173: merr = dmmoab->mbiface->tag_iterate(vtag, dmmoab->vlocal->begin(), dmmoab->vlocal->end(), count, reinterpret_cast<void *&>(marray), false);
174: MBERRNM(merr);
175: 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);
177: i = 0;
178: for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
179: for (f = 0; f < dmmoab->numFields; f++, i++) (*varray)[dmmoab->lidmap[(PetscInt)*iter - dmmoab->seqstart] * dmmoab->numFields + f] = marray[i];
180: //(*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]]=marray[i];
181: }
182: }
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /*@C
187: DMMoabVecRestoreArray - Restores the writable direct access array obtained via DMMoabVecGetArray
189: Collective
191: Input Parameters:
192: + dm - The DMMoab object being set
193: . vec - The Vector whose underlying data is requested
194: - array - The local data array
196: Level: intermediate
198: .seealso: `DMMoabVecGetArray()`, `DMMoabVecGetArrayRead()`, `DMMoabVecRestoreArrayRead()`
199: @*/
200: PetscErrorCode DMMoabVecRestoreArray(DM dm, Vec vec, void *array)
201: {
202: DM_Moab *dmmoab;
203: moab::ErrorCode merr;
204: moab::Tag vtag;
205: PetscInt count, i, f;
206: PetscScalar **varray;
207: PetscScalar *marray;
208: PetscContainer moabdata;
209: Vec_MOAB *vmoab, *xmoab;
211: PetscFunctionBegin;
214: PetscAssertPointer(array, 3);
215: dmmoab = (DM_Moab *)dm->data;
217: /* Get the Vec_MOAB struct for the original vector */
218: PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
219: PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
221: /* Get the real scalar array handle */
222: varray = reinterpret_cast<PetscScalar **>(array);
224: if (vmoab->is_native_vec) {
225: /* Get the Vec_MOAB struct for the original vector */
226: PetscCall(PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject *)&moabdata));
227: PetscCall(PetscContainerGetPointer(moabdata, &xmoab));
229: /* get the local representation of the arrays from Vectors */
230: PetscCall(VecRestoreArray(xmoab->local, varray));
231: PetscCall(VecGhostUpdateBegin(vmoab->local, ADD_VALUES, SCATTER_REVERSE));
232: PetscCall(VecGhostUpdateEnd(vmoab->local, ADD_VALUES, SCATTER_REVERSE));
233: PetscCall(VecGhostRestoreLocalForm(vmoab->local, &xmoab->local));
235: /* restore local pieces */
236: PetscCall(DMLocalToGlobalBegin(dm, vmoab->local, INSERT_VALUES, vec));
237: PetscCall(DMLocalToGlobalEnd(dm, vmoab->local, INSERT_VALUES, vec));
238: PetscCall(DMRestoreLocalVector(dm, &vmoab->local));
239: } else {
240: /* Get the MOAB private data */
241: PetscCall(DMMoabGetVecTag(vec, &vtag));
243: /* Get the array data for local entities */
244: merr = dmmoab->mbiface->tag_iterate(vtag, dmmoab->vlocal->begin(), dmmoab->vlocal->end(), count, reinterpret_cast<void *&>(marray), false);
245: MBERRNM(merr);
246: 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);
248: i = 0;
249: for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
250: for (f = 0; f < dmmoab->numFields; f++, i++) marray[i] = (*varray)[dmmoab->lidmap[(PetscInt)*iter - dmmoab->seqstart] * dmmoab->numFields + f];
251: //marray[i] = (*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]];
252: }
254: #ifdef MOAB_HAVE_MPI
255: /* reduce the tags correctly -> should probably let the user choose how to reduce in the future
256: For all FEM residual based assembly calculations, MPI_SUM should serve well */
257: merr = dmmoab->pcomm->reduce_tags(vtag, MPI_SUM, *dmmoab->vlocal);
258: MBERRV(dmmoab->mbiface, merr);
259: #endif
260: PetscCall(PetscFree(*varray));
261: }
262: PetscFunctionReturn(PETSC_SUCCESS);
263: }
265: /*@C
266: 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
268: Collective
270: Input Parameters:
271: + dm - The DMMoab object being set
272: - vec - The Vector whose underlying data is requested
274: Output Parameter:
275: . array - The local data array
277: Level: intermediate
279: .seealso: `DMMoabVecRestoreArrayRead()`, `DMMoabVecGetArray()`, `DMMoabVecRestoreArray()`
280: @*/
281: PetscErrorCode DMMoabVecGetArrayRead(DM dm, Vec vec, void *array)
282: {
283: DM_Moab *dmmoab;
284: moab::ErrorCode merr;
285: PetscInt count, i, f;
286: moab::Tag vtag;
287: PetscScalar **varray;
288: PetscScalar *marray;
289: PetscContainer moabdata;
290: Vec_MOAB *vmoab, *xmoab;
292: PetscFunctionBegin;
295: PetscAssertPointer(array, 3);
296: dmmoab = (DM_Moab *)dm->data;
298: /* Get the Vec_MOAB struct for the original vector */
299: PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
300: PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
302: /* Get the real scalar array handle */
303: varray = reinterpret_cast<PetscScalar **>(array);
305: if (vmoab->is_native_vec) {
306: /* get the local representation of the arrays from Vectors */
307: PetscCall(DMGetLocalVector(dm, &vmoab->local));
308: PetscCall(DMGlobalToLocalBegin(dm, vec, INSERT_VALUES, vmoab->local));
309: PetscCall(DMGlobalToLocalEnd(dm, vec, INSERT_VALUES, vmoab->local));
311: /* Get the Vec_MOAB struct for the original vector */
312: PetscCall(PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject *)&moabdata));
313: PetscCall(PetscContainerGetPointer(moabdata, &xmoab));
315: /* get the local representation of the arrays from Vectors */
316: PetscCall(VecGhostGetLocalForm(vmoab->local, &xmoab->local));
317: PetscCall(VecGhostUpdateBegin(vmoab->local, INSERT_VALUES, SCATTER_FORWARD));
318: PetscCall(VecGhostUpdateEnd(vmoab->local, INSERT_VALUES, SCATTER_FORWARD));
319: PetscCall(VecGetArray(xmoab->local, varray));
320: } else {
321: /* Get the MOAB private data */
322: PetscCall(DMMoabGetVecTag(vec, &vtag));
324: #ifdef MOAB_HAVE_MPI
325: /* exchange the data into ghost cells first */
326: merr = dmmoab->pcomm->exchange_tags(vtag, *dmmoab->vlocal);
327: MBERRNM(merr);
328: #endif
329: PetscCall(PetscMalloc1((dmmoab->nloc + dmmoab->nghost) * dmmoab->numFields, varray));
331: /* Get the array data for local entities */
332: merr = dmmoab->mbiface->tag_iterate(vtag, dmmoab->vlocal->begin(), dmmoab->vlocal->end(), count, reinterpret_cast<void *&>(marray), false);
333: MBERRNM(merr);
334: 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);
336: i = 0;
337: for (moab::Range::iterator iter = dmmoab->vlocal->begin(); iter != dmmoab->vlocal->end(); iter++) {
338: for (f = 0; f < dmmoab->numFields; f++, i++) (*varray)[dmmoab->lidmap[(PetscInt)*iter - dmmoab->seqstart] * dmmoab->numFields + f] = marray[i];
339: //(*varray)[dmmoab->llmap[dmmoab->lidmap[((PetscInt)*iter-dmmoab->seqstart)]*dmmoab->numFields+f]]=marray[i];
340: }
341: }
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
345: /*@C
346: DMMoabVecRestoreArrayRead - Restores the read-only direct access array obtained via DMMoabVecGetArray
348: Collective
350: Input Parameters:
351: + dm - The DMMoab object being set
352: . vec - The Vector whose underlying data is requested
353: - array - The local data array
355: Level: intermediate
357: .seealso: `DMMoabVecGetArrayRead()`, `DMMoabVecGetArray()`, `DMMoabVecRestoreArray()`
358: @*/
359: PetscErrorCode DMMoabVecRestoreArrayRead(DM dm, Vec vec, void *array)
360: {
361: PetscScalar **varray;
362: PetscContainer moabdata;
363: Vec_MOAB *vmoab, *xmoab;
365: PetscFunctionBegin;
368: PetscAssertPointer(array, 3);
370: /* Get the Vec_MOAB struct for the original vector */
371: PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
372: PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
374: /* Get the real scalar array handle */
375: varray = reinterpret_cast<PetscScalar **>(array);
377: if (vmoab->is_native_vec) {
378: /* Get the Vec_MOAB struct for the original vector */
379: PetscCall(PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject *)&moabdata));
380: PetscCall(PetscContainerGetPointer(moabdata, &xmoab));
382: /* restore the local representation of the arrays from Vectors */
383: PetscCall(VecRestoreArray(xmoab->local, varray));
384: PetscCall(VecGhostRestoreLocalForm(vmoab->local, &xmoab->local));
386: /* restore local pieces */
387: PetscCall(DMRestoreLocalVector(dm, &vmoab->local));
388: } else {
389: /* Nothing to do but just free the memory allocated before */
390: PetscCall(PetscFree(*varray));
391: }
392: PetscFunctionReturn(PETSC_SUCCESS);
393: }
395: static PetscErrorCode DMCreateVector_Moab_Private(DM dm, moab::Tag tag, const moab::Range *userrange, PetscBool is_global_vec, PetscBool destroy_tag, Vec *vec)
396: {
397: moab::ErrorCode merr;
398: PetscBool is_newtag;
399: const moab::Range *range;
400: PetscInt count, lnative_vec, gnative_vec;
401: std::string ttname;
402: PetscScalar *data_ptr, *defaultvals;
404: Vec_MOAB *vmoab;
405: DM_Moab *dmmoab = (DM_Moab *)dm->data;
406: #ifdef MOAB_HAVE_MPI
407: moab::ParallelComm *pcomm = dmmoab->pcomm;
408: #endif
409: moab::Interface *mbiface = dmmoab->mbiface;
411: PetscFunctionBegin;
412: PetscCheck(sizeof(PetscReal) == sizeof(PetscScalar), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "MOAB tags only support Real types (Complex-type unsupported)");
413: if (!userrange) range = dmmoab->vowned;
414: else range = userrange;
415: PetscCheck(range, PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Input range cannot be empty or call DMSetUp first.");
417: #if !defined(USE_NATIVE_PETSCVEC)
418: /* If the tag data is in a single sequence, use PETSc native vector since tag_iterate isn't useful anymore */
419: lnative_vec = (range->psize() - 1);
420: #else
421: lnative_vec = 1; /* NOTE: Testing PETSc vector: will force to create native vector all the time */
422: // lnative_vec=0; /* NOTE: Testing MOAB vector: will force to create MOAB tag_iterate based vector all the time */
423: #endif
425: #ifdef MOAB_HAVE_MPI
426: PetscCallMPI(MPIU_Allreduce(&lnative_vec, &gnative_vec, 1, MPI_INT, MPI_MAX, ((PetscObject)dm)->comm));
427: #else
428: gnative_vec = lnative_vec;
429: #endif
431: /* Create the MOAB internal data object */
432: PetscCall(PetscNew(&vmoab));
433: vmoab->is_native_vec = (gnative_vec > 0 ? PETSC_TRUE : PETSC_FALSE);
435: if (!vmoab->is_native_vec) {
436: merr = moab::MB_SUCCESS;
437: if (tag != 0) merr = mbiface->tag_get_name(tag, ttname);
438: if (!ttname.length() || merr != moab::MB_SUCCESS) {
439: /* get the new name for the anonymous MOABVec -> the tag_name will be destroyed along with Tag */
440: char *tag_name = NULL;
441: #ifdef MOAB_HAVE_MPI
442: PetscCall(DMVecCreateTagName_Moab_Private(mbiface, pcomm, &tag_name));
443: #else
444: PetscCall(DMVecCreateTagName_Moab_Private(mbiface, &tag_name));
445: #endif
446: is_newtag = PETSC_TRUE;
448: /* Create the default value for the tag (all zeros) */
449: PetscCall(PetscCalloc1(dmmoab->numFields, &defaultvals));
451: /* Create the tag */
452: merr = mbiface->tag_get_handle(tag_name, dmmoab->numFields, moab::MB_TYPE_DOUBLE, tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT, defaultvals);
453: MBERRNM(merr);
454: PetscCall(PetscFree(tag_name));
455: PetscCall(PetscFree(defaultvals));
456: } else {
457: /* Make sure the tag data is of type "double" */
458: moab::DataType tag_type;
459: merr = mbiface->tag_get_data_type(tag, tag_type);
460: MBERRNM(merr);
461: PetscCheck(tag_type == moab::MB_TYPE_DOUBLE, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Tag data type must be MB_TYPE_DOUBLE");
462: is_newtag = destroy_tag;
463: }
465: vmoab->tag = tag;
466: vmoab->new_tag = is_newtag;
467: }
468: vmoab->mbiface = mbiface;
469: #ifdef MOAB_HAVE_MPI
470: vmoab->pcomm = pcomm;
471: #endif
472: vmoab->is_global_vec = is_global_vec;
473: vmoab->tag_size = dmmoab->bs;
475: if (vmoab->is_native_vec) {
476: /* Create the PETSc Vector directly and attach our functions accordingly */
477: if (!is_global_vec) {
478: /* This is an MPI Vector with ghosted padding */
479: PetscCall(VecCreateGhostBlock(((PetscObject)dm)->comm, dmmoab->bs, dmmoab->numFields * dmmoab->nloc, dmmoab->numFields * dmmoab->n, dmmoab->nghost, &dmmoab->gsindices[dmmoab->nloc], vec));
480: } else {
481: /* This is an MPI/SEQ Vector */
482: PetscCall(VecCreate(((PetscObject)dm)->comm, vec));
483: PetscCall(VecSetSizes(*vec, dmmoab->numFields * dmmoab->nloc, PETSC_DECIDE));
484: PetscCall(VecSetBlockSize(*vec, dmmoab->bs));
485: PetscCall(VecSetType(*vec, VECMPI));
486: }
487: } else {
488: /* Call tag_iterate. This will cause MOAB to allocate memory for the
489: tag data if it hasn't already happened */
490: merr = mbiface->tag_iterate(tag, range->begin(), range->end(), count, (void *&)data_ptr);
491: MBERRNM(merr);
493: /* set the reference for vector range */
494: vmoab->tag_range = new moab::Range(*range);
495: merr = mbiface->tag_get_length(tag, dmmoab->numFields);
496: MBERRNM(merr);
498: /* Create the PETSc Vector
499: Query MOAB mesh to check if there are any ghosted entities
500: -> if we do, create a ghosted vector to map correctly to the same layout
501: -> else, create a non-ghosted parallel vector */
502: if (!is_global_vec) {
503: /* This is an MPI Vector with ghosted padding */
504: PetscCall(VecCreateGhostBlockWithArray(((PetscObject)dm)->comm, dmmoab->bs, dmmoab->numFields * dmmoab->nloc, dmmoab->numFields * dmmoab->n, dmmoab->nghost, &dmmoab->gsindices[dmmoab->nloc], data_ptr, vec));
505: } else {
506: /* This is an MPI Vector without ghosted padding */
507: PetscCall(VecCreateMPIWithArray(((PetscObject)dm)->comm, dmmoab->bs, dmmoab->numFields * range->size(), PETSC_DECIDE, data_ptr, vec));
508: }
509: }
510: PetscCall(VecSetFromOptions(*vec));
512: /* create a container and store the internal MOAB data for faster access based on Entities etc */
513: PetscCall(PetscObjectContainerCompose((PetscObject)*vec, "MOABData", vmoab, DMVecCtxDestroy_Moab));
514: (*vec)->ops->duplicate = DMVecDuplicate_Moab;
516: /* Vector created, manually set local to global mapping */
517: if (dmmoab->ltog_map) PetscCall(VecSetLocalToGlobalMapping(*vec, dmmoab->ltog_map));
519: /* set the DM reference to the vector */
520: PetscCall(VecSetDM(*vec, dm));
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: /* DMVecCreateTagName_Moab_Private
525: *
526: * Creates a unique tag name that will be shared across processes. If
527: * pcomm is NULL, then this is a serial vector. A unique tag name
528: * will be returned in tag_name in either case.
529: *
530: * The tag names have the format _PETSC_VEC_N where N is some integer.
531: *
532: * NOTE: The tag_name is allocated in this routine; The user needs to free
533: * the character array.
534: */
535: #ifdef MOAB_HAVE_MPI
536: PetscErrorCode DMVecCreateTagName_Moab_Private(moab::Interface *mbiface, moab::ParallelComm *pcomm, char **tag_name)
537: #else
538: static PetscErrorCode DMVecCreateTagName_Moab_Private(moab::Interface *mbiface, char **tag_name)
539: #endif
540: {
541: moab::ErrorCode mberr;
542: PetscInt n, global_n;
543: moab::Tag indexTag;
545: PetscFunctionBegin;
546: const char *PVEC_PREFIX = "__PETSC_VEC_";
547: PetscCall(PetscMalloc1(PETSC_MAX_PATH_LEN, tag_name));
549: moab::EntityHandle rootset = mbiface->get_root_set();
551: /* Check to see if there are any PETSc vectors defined */
552: /* Create a tag in MOAB mesh to index and keep track of number of PETSc vec tags */
553: mberr = mbiface->tag_get_handle("__PETSC_VECS__", 1, moab::MB_TYPE_INTEGER, indexTag, moab::MB_TAG_SPARSE | moab::MB_TAG_CREAT, 0);
554: MBERRNM(mberr);
555: mberr = mbiface->tag_get_data(indexTag, &rootset, 1, &n);
556: if (mberr == moab::MB_TAG_NOT_FOUND) n = 0; /* this is the first temporary vector */
557: else MBERRNM(mberr);
559: /* increment the new value of n */
560: ++n;
562: #ifdef MOAB_HAVE_MPI
563: /* Make sure that n is consistent across all processes */
564: PetscCallMPI(MPIU_Allreduce(&n, &global_n, 1, MPI_INT, MPI_MAX, pcomm->comm()));
565: #else
566: global_n = n;
567: #endif
569: /* Set the new name accordingly and return */
570: PetscCall(PetscSNPrintf(*tag_name, PETSC_MAX_PATH_LEN - 1, "%s_%" PetscInt_FMT, PVEC_PREFIX, global_n));
571: mberr = mbiface->tag_set_data(indexTag, &rootset, 1, (const void *)&global_n);
572: MBERRNM(mberr);
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM dm, Vec *gvec)
577: {
578: DM_Moab *dmmoab = (DM_Moab *)dm->data;
580: PetscFunctionBegin;
582: PetscAssertPointer(gvec, 2);
583: PetscCall(DMCreateVector_Moab_Private(dm, NULL, dmmoab->vowned, PETSC_TRUE, PETSC_TRUE, gvec));
584: PetscFunctionReturn(PETSC_SUCCESS);
585: }
587: PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM dm, Vec *lvec)
588: {
589: DM_Moab *dmmoab = (DM_Moab *)dm->data;
591: PetscFunctionBegin;
593: PetscAssertPointer(lvec, 2);
594: PetscCall(DMCreateVector_Moab_Private(dm, NULL, dmmoab->vlocal, PETSC_FALSE, PETSC_TRUE, lvec));
595: PetscFunctionReturn(PETSC_SUCCESS);
596: }
598: static PetscErrorCode DMVecDuplicate_Moab(Vec x, Vec *y)
599: {
600: DM dm;
601: PetscContainer moabdata;
602: Vec_MOAB *vmoab;
604: PetscFunctionBegin;
606: PetscAssertPointer(y, 2);
608: /* Get the Vec_MOAB struct for the original vector */
609: PetscCall(PetscObjectQuery((PetscObject)x, "MOABData", (PetscObject *)&moabdata));
610: PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
612: PetscCall(VecGetDM(x, &dm));
615: PetscCall(DMCreateVector_Moab_Private(dm, NULL, vmoab->tag_range, vmoab->is_global_vec, PETSC_TRUE, y));
616: PetscCall(VecSetDM(*y, dm));
617: PetscFunctionReturn(PETSC_SUCCESS);
618: }
620: static PetscErrorCode DMVecCtxDestroy_Moab(PetscCtxRt ctx)
621: {
622: Vec_MOAB *vmoab = *(Vec_MOAB **)ctx;
623: moab::ErrorCode merr;
625: PetscFunctionBegin;
626: if (vmoab->new_tag && vmoab->tag) {
627: /* Tag was created via a call to VecDuplicate, delete the underlying tag in MOAB */
628: merr = vmoab->mbiface->tag_delete(vmoab->tag);
629: MBERRNM(merr);
630: }
631: delete vmoab->tag_range;
632: vmoab->tag = NULL;
633: vmoab->mbiface = NULL;
634: #ifdef MOAB_HAVE_MPI
635: vmoab->pcomm = NULL;
636: #endif
637: PetscCall(PetscFree(vmoab));
638: PetscFunctionReturn(PETSC_SUCCESS);
639: }
641: PETSC_EXTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM dm, Vec g, InsertMode mode, Vec l)
642: {
643: DM_Moab *dmmoab = (DM_Moab *)dm->data;
645: PetscFunctionBegin;
646: PetscCall(VecScatterBegin(dmmoab->ltog_sendrecv, g, l, mode, SCATTER_REVERSE));
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: PETSC_EXTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM dm, Vec g, InsertMode mode, Vec l)
651: {
652: DM_Moab *dmmoab = (DM_Moab *)dm->data;
654: PetscFunctionBegin;
655: PetscCall(VecScatterEnd(dmmoab->ltog_sendrecv, g, l, mode, SCATTER_REVERSE));
656: PetscFunctionReturn(PETSC_SUCCESS);
657: }
659: PETSC_EXTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM dm, Vec l, InsertMode mode, Vec g)
660: {
661: DM_Moab *dmmoab = (DM_Moab *)dm->data;
663: PetscFunctionBegin;
664: PetscCall(VecScatterBegin(dmmoab->ltog_sendrecv, l, g, mode, SCATTER_FORWARD));
665: PetscFunctionReturn(PETSC_SUCCESS);
666: }
668: PETSC_EXTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM dm, Vec l, InsertMode mode, Vec g)
669: {
670: DM_Moab *dmmoab = (DM_Moab *)dm->data;
672: PetscFunctionBegin;
673: PetscCall(VecScatterEnd(dmmoab->ltog_sendrecv, l, g, mode, SCATTER_FORWARD));
674: PetscFunctionReturn(PETSC_SUCCESS);
675: }