1 #include <petsc/private/dmmbimpl.h> /*I "petscdmmoab.h" I*/
2 #include <petsc/private/vecimpl.h>
3
4 #include <petscdmmoab.h>
5 #include <MBTagConventions.hpp>
6
7 #define USE_NATIVE_PETSCVEC
8
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
18
19 /*@C
20 DMMoabCreateVector - Create a Vec from either an existing tag, or a specified tag size, and a range of entities
21
22 Collective
23
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
30
31 Output Parameter:
32 . vec - The created vector
33
34 Level: beginner
35
36 .seealso: `VecCreate()`
37 @*/
DMMoabCreateVector(DM dm,moab::Tag tag,const moab::Range * range,PetscBool is_global_vec,PetscBool destroy_tag,Vec * vec)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.");
42
43 PetscCall(DMCreateVector_Moab_Private(dm, tag, range, is_global_vec, destroy_tag, vec));
44 PetscFunctionReturn(PETSC_SUCCESS);
45 }
46
47 /*@C
48 DMMoabGetVecTag - Get the MOAB tag associated with this Vec
49
50 Input Parameter:
51 . vec - Vec being queried
52
53 Output Parameter:
54 . tag - Tag associated with this Vec. NULL if vec is a native PETSc Vec.
55
56 Level: beginner
57
58 .seealso: `DMMoabCreateVector()`, `DMMoabGetVecRange()`
59 @*/
DMMoabGetVecTag(Vec vec,moab::Tag * tag)60 PetscErrorCode DMMoabGetVecTag(Vec vec, moab::Tag *tag)
61 {
62 PetscContainer moabdata;
63 Vec_MOAB *vmoab;
64
65 PetscFunctionBegin;
66 PetscAssertPointer(tag, 2);
67
68 /* Get the MOAB private data */
69 PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
70 PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
71
72 *tag = vmoab->tag;
73 PetscFunctionReturn(PETSC_SUCCESS);
74 }
75
76 /*@C
77 DMMoabGetVecRange - Get the MOAB entities associated with this Vec
78
79 Input Parameter:
80 . vec - Vec being queried
81
82 Output Parameter:
83 . range - Entities associated with this Vec. NULL if vec is a native PETSc Vec.
84
85 Level: beginner
86
87 .seealso: `DMMoabCreateVector()`, `DMMoabGetVecTag()`
88 @*/
DMMoabGetVecRange(Vec vec,moab::Range * range)89 PetscErrorCode DMMoabGetVecRange(Vec vec, moab::Range *range)
90 {
91 PetscContainer moabdata;
92 Vec_MOAB *vmoab;
93
94 PetscFunctionBegin;
95 PetscAssertPointer(range, 2);
96
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 }
103
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
106
107 Collective
108
109 Input Parameters:
110 + dm - The DMMoab object being set
111 - vec - The Vector whose underlying data is requested
112
113 Output Parameter:
114 . array - The local data array
115
116 Level: intermediate
117
118 .seealso: `DMMoabVecRestoreArray()`, `DMMoabVecGetArrayRead()`, `DMMoabVecRestoreArrayRead()`
119 @*/
DMMoabVecGetArray(DM dm,Vec vec,void * array)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;
130
131 PetscFunctionBegin;
132 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
133 PetscValidHeaderSpecific(vec, VEC_CLASSID, 2);
134 PetscAssertPointer(array, 3);
135 dmmoab = (DM_Moab *)dm->data;
136
137 /* Get the Vec_MOAB struct for the original vector */
138 PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
139 PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
140
141 /* Get the real scalar array handle */
142 varray = reinterpret_cast<PetscScalar **>(array);
143
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));
149
150 /* Get the Vec_MOAB struct for the original vector */
151 PetscCall(PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject *)&moabdata));
152 PetscCall(PetscContainerGetPointer(moabdata, &xmoab));
153
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));
158
159 PetscCall(VecGetArray(xmoab->local, varray));
160 } else {
161 /* Get the MOAB private data */
162 PetscCall(DMMoabGetVecTag(vec, &vtag));
163
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
169
170 PetscCall(PetscMalloc1((dmmoab->nloc + dmmoab->nghost) * dmmoab->numFields, varray));
171
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);
176
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 }
185
186 /*@C
187 DMMoabVecRestoreArray - Restores the writable direct access array obtained via DMMoabVecGetArray
188
189 Collective
190
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
195
196 Level: intermediate
197
198 .seealso: `DMMoabVecGetArray()`, `DMMoabVecGetArrayRead()`, `DMMoabVecRestoreArrayRead()`
199 @*/
DMMoabVecRestoreArray(DM dm,Vec vec,void * array)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;
210
211 PetscFunctionBegin;
212 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
213 PetscValidHeaderSpecific(vec, VEC_CLASSID, 2);
214 PetscAssertPointer(array, 3);
215 dmmoab = (DM_Moab *)dm->data;
216
217 /* Get the Vec_MOAB struct for the original vector */
218 PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
219 PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
220
221 /* Get the real scalar array handle */
222 varray = reinterpret_cast<PetscScalar **>(array);
223
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));
228
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));
234
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));
242
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);
247
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 }
253
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 }
264
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
267
268 Collective
269
270 Input Parameters:
271 + dm - The DMMoab object being set
272 - vec - The Vector whose underlying data is requested
273
274 Output Parameter:
275 . array - The local data array
276
277 Level: intermediate
278
279 .seealso: `DMMoabVecRestoreArrayRead()`, `DMMoabVecGetArray()`, `DMMoabVecRestoreArray()`
280 @*/
DMMoabVecGetArrayRead(DM dm,Vec vec,void * array)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;
291
292 PetscFunctionBegin;
293 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
294 PetscValidHeaderSpecific(vec, VEC_CLASSID, 2);
295 PetscAssertPointer(array, 3);
296 dmmoab = (DM_Moab *)dm->data;
297
298 /* Get the Vec_MOAB struct for the original vector */
299 PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
300 PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
301
302 /* Get the real scalar array handle */
303 varray = reinterpret_cast<PetscScalar **>(array);
304
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));
310
311 /* Get the Vec_MOAB struct for the original vector */
312 PetscCall(PetscObjectQuery((PetscObject)vmoab->local, "MOABData", (PetscObject *)&moabdata));
313 PetscCall(PetscContainerGetPointer(moabdata, &xmoab));
314
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));
323
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));
330
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);
335
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 }
344
345 /*@C
346 DMMoabVecRestoreArrayRead - Restores the read-only direct access array obtained via DMMoabVecGetArray
347
348 Collective
349
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
354
355 Level: intermediate
356
357 .seealso: `DMMoabVecGetArrayRead()`, `DMMoabVecGetArray()`, `DMMoabVecRestoreArray()`
358 @*/
DMMoabVecRestoreArrayRead(DM dm,Vec vec,void * array)359 PetscErrorCode DMMoabVecRestoreArrayRead(DM dm, Vec vec, void *array)
360 {
361 PetscScalar **varray;
362 PetscContainer moabdata;
363 Vec_MOAB *vmoab, *xmoab;
364
365 PetscFunctionBegin;
366 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
367 PetscValidHeaderSpecific(vec, VEC_CLASSID, 2);
368 PetscAssertPointer(array, 3);
369
370 /* Get the Vec_MOAB struct for the original vector */
371 PetscCall(PetscObjectQuery((PetscObject)vec, "MOABData", (PetscObject *)&moabdata));
372 PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
373
374 /* Get the real scalar array handle */
375 varray = reinterpret_cast<PetscScalar **>(array);
376
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));
381
382 /* restore the local representation of the arrays from Vectors */
383 PetscCall(VecRestoreArray(xmoab->local, varray));
384 PetscCall(VecGhostRestoreLocalForm(vmoab->local, &xmoab->local));
385
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 }
394
DMCreateVector_Moab_Private(DM dm,moab::Tag tag,const moab::Range * userrange,PetscBool is_global_vec,PetscBool destroy_tag,Vec * vec)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;
403
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;
410
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.");
416
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
424
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
430
431 /* Create the MOAB internal data object */
432 PetscCall(PetscNew(&vmoab));
433 vmoab->is_native_vec = (gnative_vec > 0 ? PETSC_TRUE : PETSC_FALSE);
434
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;
447
448 /* Create the default value for the tag (all zeros) */
449 PetscCall(PetscCalloc1(dmmoab->numFields, &defaultvals));
450
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 }
464
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;
474
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);
492
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);
497
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));
511
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;
515
516 /* Vector created, manually set local to global mapping */
517 if (dmmoab->ltog_map) PetscCall(VecSetLocalToGlobalMapping(*vec, dmmoab->ltog_map));
518
519 /* set the DM reference to the vector */
520 PetscCall(VecSetDM(*vec, dm));
521 PetscFunctionReturn(PETSC_SUCCESS);
522 }
523
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
DMVecCreateTagName_Moab_Private(moab::Interface * mbiface,moab::ParallelComm * pcomm,char ** tag_name)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;
544
545 PetscFunctionBegin;
546 const char *PVEC_PREFIX = "__PETSC_VEC_";
547 PetscCall(PetscMalloc1(PETSC_MAX_PATH_LEN, tag_name));
548
549 moab::EntityHandle rootset = mbiface->get_root_set();
550
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);
558
559 /* increment the new value of n */
560 ++n;
561
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
568
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 }
575
DMCreateGlobalVector_Moab(DM dm,Vec * gvec)576 PETSC_EXTERN PetscErrorCode DMCreateGlobalVector_Moab(DM dm, Vec *gvec)
577 {
578 DM_Moab *dmmoab = (DM_Moab *)dm->data;
579
580 PetscFunctionBegin;
581 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
582 PetscAssertPointer(gvec, 2);
583 PetscCall(DMCreateVector_Moab_Private(dm, NULL, dmmoab->vowned, PETSC_TRUE, PETSC_TRUE, gvec));
584 PetscFunctionReturn(PETSC_SUCCESS);
585 }
586
DMCreateLocalVector_Moab(DM dm,Vec * lvec)587 PETSC_EXTERN PetscErrorCode DMCreateLocalVector_Moab(DM dm, Vec *lvec)
588 {
589 DM_Moab *dmmoab = (DM_Moab *)dm->data;
590
591 PetscFunctionBegin;
592 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
593 PetscAssertPointer(lvec, 2);
594 PetscCall(DMCreateVector_Moab_Private(dm, NULL, dmmoab->vlocal, PETSC_FALSE, PETSC_TRUE, lvec));
595 PetscFunctionReturn(PETSC_SUCCESS);
596 }
597
DMVecDuplicate_Moab(Vec x,Vec * y)598 static PetscErrorCode DMVecDuplicate_Moab(Vec x, Vec *y)
599 {
600 DM dm;
601 PetscContainer moabdata;
602 Vec_MOAB *vmoab;
603
604 PetscFunctionBegin;
605 PetscValidHeaderSpecific(x, VEC_CLASSID, 1);
606 PetscAssertPointer(y, 2);
607
608 /* Get the Vec_MOAB struct for the original vector */
609 PetscCall(PetscObjectQuery((PetscObject)x, "MOABData", (PetscObject *)&moabdata));
610 PetscCall(PetscContainerGetPointer(moabdata, &vmoab));
611
612 PetscCall(VecGetDM(x, &dm));
613 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
614
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 }
619
DMVecCtxDestroy_Moab(PetscCtxRt ctx)620 static PetscErrorCode DMVecCtxDestroy_Moab(PetscCtxRt ctx)
621 {
622 Vec_MOAB *vmoab = *(Vec_MOAB **)ctx;
623 moab::ErrorCode merr;
624
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 }
640
DMGlobalToLocalBegin_Moab(DM dm,Vec g,InsertMode mode,Vec l)641 PETSC_EXTERN PetscErrorCode DMGlobalToLocalBegin_Moab(DM dm, Vec g, InsertMode mode, Vec l)
642 {
643 DM_Moab *dmmoab = (DM_Moab *)dm->data;
644
645 PetscFunctionBegin;
646 PetscCall(VecScatterBegin(dmmoab->ltog_sendrecv, g, l, mode, SCATTER_REVERSE));
647 PetscFunctionReturn(PETSC_SUCCESS);
648 }
649
DMGlobalToLocalEnd_Moab(DM dm,Vec g,InsertMode mode,Vec l)650 PETSC_EXTERN PetscErrorCode DMGlobalToLocalEnd_Moab(DM dm, Vec g, InsertMode mode, Vec l)
651 {
652 DM_Moab *dmmoab = (DM_Moab *)dm->data;
653
654 PetscFunctionBegin;
655 PetscCall(VecScatterEnd(dmmoab->ltog_sendrecv, g, l, mode, SCATTER_REVERSE));
656 PetscFunctionReturn(PETSC_SUCCESS);
657 }
658
DMLocalToGlobalBegin_Moab(DM dm,Vec l,InsertMode mode,Vec g)659 PETSC_EXTERN PetscErrorCode DMLocalToGlobalBegin_Moab(DM dm, Vec l, InsertMode mode, Vec g)
660 {
661 DM_Moab *dmmoab = (DM_Moab *)dm->data;
662
663 PetscFunctionBegin;
664 PetscCall(VecScatterBegin(dmmoab->ltog_sendrecv, l, g, mode, SCATTER_FORWARD));
665 PetscFunctionReturn(PETSC_SUCCESS);
666 }
667
DMLocalToGlobalEnd_Moab(DM dm,Vec l,InsertMode mode,Vec g)668 PETSC_EXTERN PetscErrorCode DMLocalToGlobalEnd_Moab(DM dm, Vec l, InsertMode mode, Vec g)
669 {
670 DM_Moab *dmmoab = (DM_Moab *)dm->data;
671
672 PetscFunctionBegin;
673 PetscCall(VecScatterEnd(dmmoab->ltog_sendrecv, l, g, mode, SCATTER_FORWARD));
674 PetscFunctionReturn(PETSC_SUCCESS);
675 }
676