Actual source code: dmmbfield.cxx
1: #include <petsc/private/dmmbimpl.h>
3: #include <petscdmmoab.h>
5: /*@C
6: DMMoabSetFieldVector - Sets the vector reference that represents the solution associated
7: with a particular field component.
9: Not Collective
11: Input Parameters:
12: + dm - the discretization manager object
13: . ifield - the index of the field as set before via DMMoabSetFieldName.
14: - fvec - the Vector solution corresponding to the field (component)
16: Level: intermediate
18: .seealso: `DMMoabGetFieldName()`, `DMMoabSetGlobalFieldVector()`
19: @*/
20: PetscErrorCode DMMoabSetFieldVector(DM dm, PetscInt ifield, Vec fvec)
21: {
22: DM_Moab *dmmoab;
23: moab::Tag vtag, ntag;
24: const PetscScalar *varray;
25: PetscScalar *farray;
26: moab::ErrorCode merr;
27: std::string tag_name;
29: PetscFunctionBegin;
31: dmmoab = (DM_Moab *)(dm)->data;
33: PetscCheck(!(ifield < 0) && !(ifield >= dmmoab->numFields), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The field %d should be positive and less than %d.", ifield, dmmoab->numFields);
35: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
36: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
37: MBERRNM(merr);
39: PetscCall(DMMoabGetVecTag(fvec, &vtag));
41: merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
42: if (!tag_name.length() && merr != moab::MB_SUCCESS) {
43: PetscCall(VecGetArrayRead(fvec, &varray));
44: /* use the entity handle and the Dof index to set the right value */
45: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)varray);
46: MBERRNM(merr);
47: PetscCall(VecRestoreArrayRead(fvec, &varray));
48: } else {
49: PetscCall(PetscMalloc1(dmmoab->nloc, &farray));
50: /* we are using a MOAB Vec - directly copy the tag data to new one */
51: merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void *)farray);
52: MBERRNM(merr);
53: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)farray);
54: MBERRNM(merr);
55: /* make sure the parallel exchange for ghosts are done appropriately */
56: PetscCall(PetscFree(farray));
57: }
58: #ifdef MOAB_HAVE_MPI
59: merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vowned);
60: MBERRNM(merr);
61: #endif
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: /*@C
66: DMMoabSetGlobalFieldVector - Sets the vector reference that represents the global solution associated
67: with all fields (components) managed by DM.
68: A useful utility when updating the DM solution after a solve, to be serialized with the mesh for
69: checkpointing purposes.
71: Not Collective
73: Input Parameters:
74: + dm - the discretization manager object
75: - fvec - the global Vector solution corresponding to all the fields managed by DM
77: Level: intermediate
79: .seealso: `DMMoabGetFieldName()`, `DMMoabSetFieldVector()`
80: @*/
81: PetscErrorCode DMMoabSetGlobalFieldVector(DM dm, Vec fvec)
82: {
83: DM_Moab *dmmoab;
84: moab::Tag vtag, ntag;
85: const PetscScalar *rarray;
86: PetscScalar *varray, *farray;
87: moab::ErrorCode merr;
88: PetscInt i, ifield;
89: std::string tag_name;
90: moab::Range::iterator iter;
92: PetscFunctionBegin;
94: dmmoab = (DM_Moab *)(dm)->data;
96: /* get the Tag corresponding to the global vector - possible that there is no tag associated.. */
97: PetscCall(DMMoabGetVecTag(fvec, &vtag));
98: merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
99: PetscCall(PetscMalloc1(dmmoab->nloc, &farray));
100: if (!tag_name.length() && merr != moab::MB_SUCCESS) {
101: /* not a MOAB vector - use VecGetSubVector to get the parts as needed */
102: PetscCall(VecGetArrayRead(fvec, &rarray));
103: for (ifield = 0; ifield < dmmoab->numFields; ++ifield) {
104: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
105: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
106: MBERRNM(merr);
108: for (i = 0; i < dmmoab->nloc; i++) farray[i] = (dmmoab->bs == 1 ? rarray[ifield * dmmoab->nloc + i] : rarray[i * dmmoab->numFields + ifield]);
110: /* use the entity handle and the Dof index to set the right value */
111: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)farray);
112: MBERRNM(merr);
113: }
114: PetscCall(VecRestoreArrayRead(fvec, &rarray));
115: } else {
116: PetscCall(PetscMalloc1(dmmoab->nloc * dmmoab->numFields, &varray));
118: /* we are using a MOAB Vec - directly copy the tag data to new one */
119: merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void *)varray);
120: MBERRNM(merr);
121: for (ifield = 0; ifield < dmmoab->numFields; ++ifield) {
122: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
123: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
124: MBERRNM(merr);
126: /* we are using a MOAB Vec - directly copy the tag data to new one */
127: for (i = 0; i < dmmoab->nloc; i++) farray[i] = (dmmoab->bs == 1 ? varray[ifield * dmmoab->nloc + i] : varray[i * dmmoab->numFields + ifield]);
129: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)farray);
130: MBERRNM(merr);
132: #ifdef MOAB_HAVE_MPI
133: /* make sure the parallel exchange for ghosts are done appropriately */
134: merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vlocal);
135: MBERRNM(merr);
136: #endif
137: }
138: PetscCall(PetscFree(varray));
139: }
140: PetscCall(PetscFree(farray));
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: /*@C
145: DMMoabSetFieldNames - Sets the number of fields and their names to be managed by the DM
147: Not Collective
149: Input Parameters:
150: + dm - the discretization manager object
151: . numFields - the total number of fields
152: - fields - the array containing the names of each field (component); Can be NULL.
154: Level: intermediate
156: .seealso: `DMMoabGetFieldName()`, `DMMoabSetFieldName()`
157: @*/
158: PetscErrorCode DMMoabSetFieldNames(DM dm, PetscInt numFields, const char *fields[])
159: {
160: PetscInt i;
161: DM_Moab *dmmoab;
163: PetscFunctionBegin;
165: dmmoab = (DM_Moab *)(dm)->data;
167: /* first deallocate the existing field structure */
168: if (dmmoab->fieldNames) {
169: for (i = 0; i < dmmoab->numFields; i++) PetscCall(PetscFree(dmmoab->fieldNames[i]));
170: PetscCall(PetscFree(dmmoab->fieldNames));
171: }
173: /* now re-allocate and assign field names */
174: dmmoab->numFields = numFields;
175: PetscCall(PetscMalloc1(numFields, &dmmoab->fieldNames));
176: if (fields) {
177: for (i = 0; i < dmmoab->numFields; i++) PetscCall(PetscStrallocpy(fields[i], (char **)&dmmoab->fieldNames[i]));
178: }
179: PetscCall(DMSetNumFields(dm, numFields));
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: /*@C
184: DMMoabGetFieldName - Gets the names of individual field components in multicomponent
185: vectors associated with a DMDA.
187: Not Collective
189: Input Parameters:
190: + dm - the discretization manager object
191: - field - field number for the DMMoab (0, 1, ... dof-1), where dof indicates the
192: number of degrees of freedom per node within the DMMoab
194: Output Parameter:
195: . fieldName - the name of the field (component)
197: Level: intermediate
199: .seealso: `DMMoabSetFieldName()`, `DMMoabSetFields()`
200: @*/
201: PetscErrorCode DMMoabGetFieldName(DM dm, PetscInt field, const char **fieldName)
202: {
203: DM_Moab *dmmoab;
205: PetscFunctionBegin;
207: dmmoab = (DM_Moab *)(dm)->data;
208: PetscCheck(!(field < 0) && !(field >= dmmoab->numFields), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "DM field %d should be in [%d, %d)", field, 0, dmmoab->numFields);
210: *fieldName = dmmoab->fieldNames[field];
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: /*@C
215: DMMoabSetFieldName - Sets the name of a field (component) managed by the DM
217: Not Collective
219: Input Parameters:
220: + dm - the discretization manager object
221: . field - the field number
222: - fieldName - the field (component) name
224: Level: intermediate
225: Notes:
226: Can only be called after DMMoabSetFields supplied with correct numFields
228: .seealso: `DMMoabGetFieldName()`, `DMMoabSetFields()`
229: @*/
230: PetscErrorCode DMMoabSetFieldName(DM dm, PetscInt field, const char *fieldName)
231: {
232: DM_Moab *dmmoab;
234: PetscFunctionBegin;
238: dmmoab = (DM_Moab *)(dm)->data;
239: PetscCheck(!(field < 0) && !(field >= dmmoab->numFields), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "DM field %d should be in [%d, %d)", field, 0, dmmoab->numFields);
241: if (dmmoab->fieldNames[field]) PetscCall(PetscFree(dmmoab->fieldNames[field]));
242: PetscCall(PetscStrallocpy(fieldName, (char **)&dmmoab->fieldNames[field]));
243: PetscFunctionReturn(PETSC_SUCCESS);
244: }
246: /*@C
247: DMMoabGetFieldDof - Gets the global degree-of-freedom of a field (component) defined on a
248: particular MOAB EntityHandle.
250: Not Collective
252: Input Parameters:
253: + dm - the discretization manager object
254: . point - the MOAB EntityHandle container which holds the field degree-of-freedom values
255: - field - the field (component) index
257: Output Parameter:
258: . dof - the global degree-of-freedom index corresponding to the field in the discrete representation (Vec, Mat)
260: Level: beginner
262: .seealso: `DMMoabGetFieldDofs()`, `DMMoabGetFieldDofsLocal()`
263: @*/
264: PetscErrorCode DMMoabGetFieldDof(DM dm, moab::EntityHandle point, PetscInt field, PetscInt *dof)
265: {
266: DM_Moab *dmmoab;
268: PetscFunctionBegin;
270: dmmoab = (DM_Moab *)(dm)->data;
272: *dof = (dmmoab->bs == 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(point) - dmmoab->seqstart] + field * dmmoab->n : dmmoab->gidmap[dmmoab->mbiface->id_from_handle(point) - dmmoab->seqstart] * dmmoab->numFields + field);
273: PetscFunctionReturn(PETSC_SUCCESS);
274: }
276: /*@C
277: DMMoabGetFieldDofs - Gets the global degree-of-freedom of a field (component) defined on an
278: array of MOAB EntityHandles.
280: Not Collective
282: Input Parameters:
283: + dm - the discretization manager object
284: . npoints - the total number of Entities in the points array
285: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
286: - field - the field (component) index
288: Output Parameter:
289: . dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
291: Level: intermediate
293: .seealso: `DMMoabGetFieldDof()`, `DMMoabGetFieldDofsLocal()`
294: @*/
295: PetscErrorCode DMMoabGetFieldDofs(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt field, PetscInt *dof)
296: {
297: PetscInt i;
298: DM_Moab *dmmoab;
300: PetscFunctionBegin;
303: dmmoab = (DM_Moab *)(dm)->data;
305: if (!dof) PetscCall(PetscMalloc1(npoints, &dof));
307: /* compute the DOF based on local blocking in the fields */
308: /* We also assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
309: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
310: for (i = 0; i < npoints; ++i)
311: dof[i] = (dmmoab->bs == 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + field * dmmoab->n : dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field);
312: PetscFunctionReturn(PETSC_SUCCESS);
313: }
315: /*@C
316: DMMoabGetFieldDofsLocal - Gets the local degrees-of-freedom of a field (component) defined on an
317: array of MOAB EntityHandles.
319: Not Collective
321: Input Parameters:
322: + dm - the discretization manager object
323: . npoints - the total number of Entities in the points array
324: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
325: - field - the field (component) index
327: Output Parameter:
328: . dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
330: Level: intermediate
332: .seealso: `DMMoabGetFieldDof()`, `DMMoabGetFieldDofs()`
333: @*/
334: PetscErrorCode DMMoabGetFieldDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt field, PetscInt *dof)
335: {
336: PetscInt i;
337: DM_Moab *dmmoab;
339: PetscFunctionBegin;
342: dmmoab = (DM_Moab *)(dm)->data;
344: if (!dof) PetscCall(PetscMalloc1(npoints, &dof));
346: /* compute the DOF based on local blocking in the fields */
347: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
348: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
349: for (i = 0; i < npoints; ++i) {
350: dof[i] = (dmmoab->bs > 1 ? dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field : dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + field * dmmoab->n);
351: }
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /*@C
356: DMMoabGetDofs - Gets the global degree-of-freedom for all fields (components) defined on an
357: array of MOAB EntityHandles.
359: Not Collective
361: Input Parameters:
362: + dm - the discretization manager object
363: . npoints - the total number of Entities in the points array
364: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
366: Output Parameter:
367: . dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
369: Level: intermediate
371: .seealso: `DMMoabGetFieldDofs()`, `DMMoabGetDofsLocal()`, `DMMoabGetDofsBlocked()`
372: @*/
373: PetscErrorCode DMMoabGetDofs(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
374: {
375: PetscInt i, field, offset;
376: DM_Moab *dmmoab;
378: PetscFunctionBegin;
381: dmmoab = (DM_Moab *)(dm)->data;
383: if (!dof) PetscCall(PetscMalloc1(dmmoab->numFields * npoints, &dof));
385: /* compute the DOF based on local blocking in the fields */
386: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
387: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
388: for (field = 0; field < dmmoab->numFields; ++field) {
389: offset = field * dmmoab->n;
390: for (i = 0; i < npoints; ++i)
391: dof[i * dmmoab->numFields + field] = (dmmoab->bs > 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field : dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + offset);
392: }
393: PetscFunctionReturn(PETSC_SUCCESS);
394: }
396: /*@C
397: DMMoabGetDofsLocal - Gets the local degree-of-freedom for all fields (components) defined on an
398: array of MOAB EntityHandles.
400: Not Collective
402: Input Parameters:
403: + dm - the discretization manager object
404: . npoints - the total number of Entities in the points array
405: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
407: Output Parameter:
408: . dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
410: Level: intermediate
412: .seealso: `DMMoabGetFieldDofs()`, `DMMoabGetDofs()`, `DMMoabGetDofsBlocked()`
413: @*/
414: PetscErrorCode DMMoabGetDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
415: {
416: PetscInt i, field, offset;
417: DM_Moab *dmmoab;
419: PetscFunctionBegin;
422: dmmoab = (DM_Moab *)(dm)->data;
424: if (!dof) PetscCall(PetscMalloc1(dmmoab->numFields * npoints, &dof));
426: /* compute the DOF based on local blocking in the fields */
427: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
428: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
429: for (field = 0; field < dmmoab->numFields; ++field) {
430: offset = field * dmmoab->n;
431: for (i = 0; i < npoints; ++i)
432: dof[i * dmmoab->numFields + field] = (dmmoab->bs > 1 ? dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field : dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + offset);
433: }
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
437: /*@C
438: DMMoabGetDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
439: array of MOAB EntityHandles. It is useful when performing Blocked(Get/Set) methods in computation
440: of element residuals and assembly of the discrete systems when all fields are co-located.
442: Not Collective
444: Input Parameters:
445: + dm - the discretization manager object
446: . npoints - the total number of Entities in the points array
447: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
449: Output Parameter:
450: . dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat)
452: Level: intermediate
454: .seealso: `DMMoabGetDofsLocal()`, `DMMoabGetDofs()`, `DMMoabGetDofsBlockedLocal()`
455: @*/
456: PetscErrorCode DMMoabGetDofsBlocked(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
457: {
458: PetscInt i;
459: DM_Moab *dmmoab;
461: PetscFunctionBegin;
464: dmmoab = (DM_Moab *)(dm)->data;
466: if (!dof) PetscCall(PetscMalloc1(npoints, &dof));
468: for (i = 0; i < npoints; ++i) dof[i] = dmmoab->gidmap[(PetscInt)points[i] - dmmoab->seqstart];
469: PetscFunctionReturn(PETSC_SUCCESS);
470: }
472: /*@C
473: DMMoabGetDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
474: array of MOAB EntityHandles. It is useful when performing local Blocked(Get/Set) methods in computation
475: of element residuals and assembly of the discrete systems when all fields are co-located.
477: Not Collective
479: Input Parameters:
480: + dm - the discretization manager object
481: . npoints - the total number of Entities in the points array
482: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
484: Output Parameter:
485: . dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat)
487: Level: intermediate
489: .seealso: `DMMoabGetDofsLocal()`, `DMMoabGetDofs()`, `DMMoabGetDofsBlockedLocal()`
490: @*/
491: PetscErrorCode DMMoabGetDofsBlockedLocal(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
492: {
493: PetscInt i;
494: DM_Moab *dmmoab;
496: PetscFunctionBegin;
499: dmmoab = (DM_Moab *)(dm)->data;
501: if (!dof) PetscCall(PetscMalloc1(npoints, &dof));
503: for (i = 0; i < npoints; ++i) dof[i] = dmmoab->lidmap[(PetscInt)points[i] - dmmoab->seqstart];
504: PetscFunctionReturn(PETSC_SUCCESS);
505: }
507: /*@C
508: DMMoabGetVertexDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
509: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
510: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
512: Not Collective
514: Input Parameters:
515: . dm - the discretization manager object
517: Output Parameter:
518: . dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
520: Level: intermediate
522: .seealso: `DMMoabGetVertexDofsBlockedLocal()`, `DMMoabGetDofsBlocked()`, `DMMoabGetDofsBlockedLocal()`
523: @*/
524: PetscErrorCode DMMoabGetVertexDofsBlocked(DM dm, PetscInt **dof)
525: {
526: DM_Moab *dmmoab;
528: PetscFunctionBegin;
530: dmmoab = (DM_Moab *)(dm)->data;
532: *dof = dmmoab->gidmap;
533: PetscFunctionReturn(PETSC_SUCCESS);
534: }
536: /*@C
537: DMMoabGetVertexDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
538: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
539: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
541: Not Collective
543: Input Parameters:
544: . dm - the discretization manager object
546: Output Parameter:
547: . dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
549: Level: intermediate
551: .seealso: `DMMoabGetVertexDofsBlocked()`, `DMMoabGetDofsBlocked()`, `DMMoabGetDofsBlockedLocal()`
552: @*/
553: PetscErrorCode DMMoabGetVertexDofsBlockedLocal(DM dm, PetscInt **dof)
554: {
555: DM_Moab *dmmoab;
557: PetscFunctionBegin;
560: dmmoab = (DM_Moab *)(dm)->data;
562: *dof = dmmoab->lidmap;
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }