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: PetscCallMOAB(dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT));
37: PetscCall(DMMoabGetVecTag(fvec, &vtag));
39: merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
40: if (!tag_name.length() && merr != moab::MB_SUCCESS) {
41: PetscCall(VecGetArrayRead(fvec, &varray));
42: /* use the entity handle and the Dof index to set the right value */
43: PetscCallMOAB(dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)varray));
44: PetscCall(VecRestoreArrayRead(fvec, &varray));
45: } else {
46: PetscCall(PetscMalloc1(dmmoab->nloc, &farray));
47: /* we are using a MOAB Vec - directly copy the tag data to new one */
48: PetscCallMOAB(dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void *)farray));
49: PetscCallMOAB(dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)farray));
50: /* make sure the parallel exchange for ghosts are done appropriately */
51: PetscCall(PetscFree(farray));
52: }
53: #ifdef MOAB_HAVE_MPI
54: PetscCallMOAB(dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vowned));
55: #endif
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: /*@C
60: DMMoabSetGlobalFieldVector - Sets the vector reference that represents the global solution associated
61: with all fields (components) managed by DM.
62: A useful utility when updating the DM solution after a solve, to be serialized with the mesh for
63: checkpointing purposes.
65: Not Collective
67: Input Parameters:
68: + dm - the discretization manager object
69: - fvec - the global Vector solution corresponding to all the fields managed by DM
71: Level: intermediate
73: .seealso: `DMMoabGetFieldName()`, `DMMoabSetFieldVector()`
74: @*/
75: PetscErrorCode DMMoabSetGlobalFieldVector(DM dm, Vec fvec)
76: {
77: DM_Moab *dmmoab;
78: moab::Tag vtag, ntag;
79: const PetscScalar *rarray;
80: PetscScalar *varray, *farray;
81: moab::ErrorCode merr;
82: PetscInt i, ifield;
83: std::string tag_name;
84: moab::Range::iterator iter;
86: PetscFunctionBegin;
88: dmmoab = (DM_Moab *)dm->data;
90: /* get the Tag corresponding to the global vector - possible that there is no tag associated.. */
91: PetscCall(DMMoabGetVecTag(fvec, &vtag));
92: merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
93: PetscCall(PetscMalloc1(dmmoab->nloc, &farray));
94: if (!tag_name.length() && merr != moab::MB_SUCCESS) {
95: /* not a MOAB vector - use VecGetSubVector to get the parts as needed */
96: PetscCall(VecGetArrayRead(fvec, &rarray));
97: for (ifield = 0; ifield < dmmoab->numFields; ++ifield) {
98: /* Create a tag in MOAB mesh to index and keep track of number of PETSc vec tags */
99: PetscCallMOAB(dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT));
100: for (i = 0; i < dmmoab->nloc; i++) farray[i] = (dmmoab->bs == 1 ? rarray[ifield * dmmoab->nloc + i] : rarray[i * dmmoab->numFields + ifield]);
102: /* use the entity handle and the Dof index to set the right value */
103: PetscCallMOAB(dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)farray));
104: }
105: PetscCall(VecRestoreArrayRead(fvec, &rarray));
106: } else {
107: PetscCall(PetscMalloc1(dmmoab->nloc * dmmoab->numFields, &varray));
109: /* we are using a MOAB Vec - directly copy the tag data to new one */
110: PetscCallMOAB(dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void *)varray));
111: for (ifield = 0; ifield < dmmoab->numFields; ++ifield) {
112: /* Create a tag in MOAB mesh to index and keep track of number of PETSc vec tags */
113: PetscCallMOAB(dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT));
115: /* we are using a MOAB Vec - directly copy the tag data to new one */
116: for (i = 0; i < dmmoab->nloc; i++) farray[i] = (dmmoab->bs == 1 ? varray[ifield * dmmoab->nloc + i] : varray[i * dmmoab->numFields + ifield]);
118: PetscCallMOAB(dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)farray));
120: #ifdef MOAB_HAVE_MPI
121: /* make sure the parallel exchange for ghosts are done appropriately */
122: PetscCallMOAB(dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vlocal));
123: #endif
124: }
125: PetscCall(PetscFree(varray));
126: }
127: PetscCall(PetscFree(farray));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: /*@C
132: DMMoabSetFieldNames - Sets the number of fields and their names to be managed by the DM
134: Not Collective
136: Input Parameters:
137: + dm - the discretization manager object
138: . numFields - the total number of fields
139: - fields - the array containing the names of each field (component); Can be `NULL`.
141: Level: intermediate
143: .seealso: `DMMoabGetFieldName()`, `DMMoabSetFieldName()`
144: @*/
145: PetscErrorCode DMMoabSetFieldNames(DM dm, PetscInt numFields, const char *fields[])
146: {
147: PetscInt i;
148: DM_Moab *dmmoab;
150: PetscFunctionBegin;
152: dmmoab = (DM_Moab *)dm->data;
154: /* first deallocate the existing field structure */
155: if (dmmoab->fieldNames) {
156: for (i = 0; i < dmmoab->numFields; i++) PetscCall(PetscFree(dmmoab->fieldNames[i]));
157: PetscCall(PetscFree(dmmoab->fieldNames));
158: }
160: /* now re-allocate and assign field names */
161: dmmoab->numFields = numFields;
162: PetscCall(PetscMalloc1(numFields, &dmmoab->fieldNames));
163: if (fields) {
164: for (i = 0; i < dmmoab->numFields; i++) PetscCall(PetscStrallocpy(fields[i], (char **)&dmmoab->fieldNames[i]));
165: }
166: PetscCall(DMSetNumFields(dm, numFields));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: /*@C
171: DMMoabGetFieldName - Gets the names of individual field components in multicomponent
172: vectors associated with a DMDA.
174: Not Collective
176: Input Parameters:
177: + dm - the discretization manager object
178: - field - field number for the DMMoab (0, 1, ... dof-1), where dof indicates the
179: number of degrees of freedom per node within the DMMoab
181: Output Parameter:
182: . fieldName - the name of the field (component)
184: Level: intermediate
186: .seealso: `DMMoabSetFieldName()`, `DMMoabSetFields()`
187: @*/
188: PetscErrorCode DMMoabGetFieldName(DM dm, PetscInt field, const char *fieldName[])
189: {
190: DM_Moab *dmmoab;
192: PetscFunctionBegin;
194: dmmoab = (DM_Moab *)dm->data;
195: 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);
197: *fieldName = dmmoab->fieldNames[field];
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: /*@C
202: DMMoabSetFieldName - Sets the name of a field (component) managed by the DM
204: Not Collective
206: Input Parameters:
207: + dm - the discretization manager object
208: . field - the field number
209: - fieldName - the field (component) name
211: Level: intermediate
213: Notes:
214: Can only be called after DMMoabSetFields supplied with correct numFields
216: .seealso: `DMMoabGetFieldName()`, `DMMoabSetFields()`
217: @*/
218: PetscErrorCode DMMoabSetFieldName(DM dm, PetscInt field, const char *fieldName)
219: {
220: DM_Moab *dmmoab;
222: PetscFunctionBegin;
224: PetscAssertPointer(fieldName, 3);
226: dmmoab = (DM_Moab *)dm->data;
227: 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);
229: PetscCall(PetscFree(dmmoab->fieldNames[field]));
230: PetscCall(PetscStrallocpy(fieldName, (char **)&dmmoab->fieldNames[field]));
231: PetscFunctionReturn(PETSC_SUCCESS);
232: }
234: /*@C
235: DMMoabGetFieldDof - Gets the global degree-of-freedom of a field (component) defined on a
236: particular MOAB EntityHandle.
238: Not Collective
240: Input Parameters:
241: + dm - the discretization manager object
242: . point - the MOAB EntityHandle container which holds the field degree-of-freedom values
243: - field - the field (component) index
245: Output Parameter:
246: . dof - the global degree-of-freedom index corresponding to the field in the discrete representation (Vec, Mat)
248: Level: beginner
250: .seealso: `DMMoabGetFieldDofs()`, `DMMoabGetFieldDofsLocal()`
251: @*/
252: PetscErrorCode DMMoabGetFieldDof(DM dm, moab::EntityHandle point, PetscInt field, PetscInt *dof)
253: {
254: DM_Moab *dmmoab;
256: PetscFunctionBegin;
258: dmmoab = (DM_Moab *)dm->data;
260: *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);
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: /*@C
265: DMMoabGetFieldDofs - Gets the global degree-of-freedom of a field (component) defined on an
266: array of MOAB EntityHandles.
268: Not Collective
270: Input Parameters:
271: + dm - the discretization manager object
272: . npoints - the total number of Entities in the points array
273: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
274: - field - the field (component) index
276: Output Parameter:
277: . dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
279: Level: intermediate
281: .seealso: `DMMoabGetFieldDof()`, `DMMoabGetFieldDofsLocal()`
282: @*/
283: PetscErrorCode DMMoabGetFieldDofs(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt field, PetscInt *dof)
284: {
285: PetscInt i;
286: DM_Moab *dmmoab;
288: PetscFunctionBegin;
290: PetscAssertPointer(points, 3);
291: dmmoab = (DM_Moab *)dm->data;
293: if (!dof) PetscCall(PetscMalloc1(npoints, &dof));
295: /* compute the DOF based on local blocking in the fields */
296: /* We also assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
297: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
298: for (i = 0; i < npoints; ++i)
299: 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);
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: /*@C
304: DMMoabGetFieldDofsLocal - Gets the local degrees-of-freedom of a field (component) defined on an
305: array of MOAB EntityHandles.
307: Not Collective
309: Input Parameters:
310: + dm - the discretization manager object
311: . npoints - the total number of Entities in the points array
312: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
313: - field - the field (component) index
315: Output Parameter:
316: . dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
318: Level: intermediate
320: .seealso: `DMMoabGetFieldDof()`, `DMMoabGetFieldDofs()`
321: @*/
322: PetscErrorCode DMMoabGetFieldDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt field, PetscInt *dof)
323: {
324: PetscInt i;
325: DM_Moab *dmmoab;
327: PetscFunctionBegin;
329: PetscAssertPointer(points, 3);
330: dmmoab = (DM_Moab *)dm->data;
332: if (!dof) PetscCall(PetscMalloc1(npoints, &dof));
334: /* compute the DOF based on local blocking in the fields */
335: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
336: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
337: for (i = 0; i < npoints; ++i) {
338: 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);
339: }
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: /*@C
344: DMMoabGetDofs - Gets the global degree-of-freedom for all fields (components) defined on an
345: array of MOAB EntityHandles.
347: Not Collective
349: Input Parameters:
350: + dm - the discretization manager object
351: . npoints - the total number of Entities in the points array
352: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
354: Output Parameter:
355: . dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
357: Level: intermediate
359: .seealso: `DMMoabGetFieldDofs()`, `DMMoabGetDofsLocal()`, `DMMoabGetDofsBlocked()`
360: @*/
361: PetscErrorCode DMMoabGetDofs(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
362: {
363: PetscInt i, field, offset;
364: DM_Moab *dmmoab;
366: PetscFunctionBegin;
368: PetscAssertPointer(points, 3);
369: dmmoab = (DM_Moab *)dm->data;
371: if (!dof) PetscCall(PetscMalloc1(dmmoab->numFields * npoints, &dof));
373: /* compute the DOF based on local blocking in the fields */
374: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
375: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
376: for (field = 0; field < dmmoab->numFields; ++field) {
377: offset = field * dmmoab->n;
378: for (i = 0; i < npoints; ++i)
379: 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);
380: }
381: PetscFunctionReturn(PETSC_SUCCESS);
382: }
384: /*@C
385: DMMoabGetDofsLocal - Gets the local degree-of-freedom for all fields (components) defined on an
386: array of MOAB EntityHandles.
388: Not Collective
390: Input Parameters:
391: + dm - the discretization manager object
392: . npoints - the total number of Entities in the points array
393: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
395: Output Parameter:
396: . dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
398: Level: intermediate
400: .seealso: `DMMoabGetFieldDofs()`, `DMMoabGetDofs()`, `DMMoabGetDofsBlocked()`
401: @*/
402: PetscErrorCode DMMoabGetDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
403: {
404: PetscInt i, field, offset;
405: DM_Moab *dmmoab;
407: PetscFunctionBegin;
409: PetscAssertPointer(points, 3);
410: dmmoab = (DM_Moab *)dm->data;
412: if (!dof) PetscCall(PetscMalloc1(dmmoab->numFields * npoints, &dof));
414: /* compute the DOF based on local blocking in the fields */
415: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
416: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
417: for (field = 0; field < dmmoab->numFields; ++field) {
418: offset = field * dmmoab->n;
419: for (i = 0; i < npoints; ++i)
420: 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);
421: }
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: /*@C
426: DMMoabGetDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
427: array of MOAB EntityHandles. It is useful when performing Blocked(Get/Set) methods in computation
428: of element residuals and assembly of the discrete systems when all fields are co-located.
430: Not Collective
432: Input Parameters:
433: + dm - the discretization manager object
434: . npoints - the total number of Entities in the points array
435: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
437: Output Parameter:
438: . dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat)
440: Level: intermediate
442: .seealso: `DMMoabGetDofsLocal()`, `DMMoabGetDofs()`, `DMMoabGetDofsBlockedLocal()`
443: @*/
444: PetscErrorCode DMMoabGetDofsBlocked(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
445: {
446: PetscInt i;
447: DM_Moab *dmmoab;
449: PetscFunctionBegin;
451: PetscAssertPointer(points, 3);
452: dmmoab = (DM_Moab *)dm->data;
454: if (!dof) PetscCall(PetscMalloc1(npoints, &dof));
456: for (i = 0; i < npoints; ++i) dof[i] = dmmoab->gidmap[(PetscInt)points[i] - dmmoab->seqstart];
457: PetscFunctionReturn(PETSC_SUCCESS);
458: }
460: /*@C
461: DMMoabGetDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
462: array of MOAB EntityHandles. It is useful when performing local Blocked(Get/Set) methods in computation
463: of element residuals and assembly of the discrete systems when all fields are co-located.
465: Not Collective
467: Input Parameters:
468: + dm - the discretization manager object
469: . npoints - the total number of Entities in the points array
470: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
472: Output Parameter:
473: . dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat)
475: Level: intermediate
477: .seealso: `DMMoabGetDofsLocal()`, `DMMoabGetDofs()`
478: @*/
479: PetscErrorCode DMMoabGetDofsBlockedLocal(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
480: {
481: PetscInt i;
482: DM_Moab *dmmoab;
484: PetscFunctionBegin;
486: PetscAssertPointer(points, 3);
487: dmmoab = (DM_Moab *)dm->data;
489: if (!dof) PetscCall(PetscMalloc1(npoints, &dof));
491: for (i = 0; i < npoints; ++i) dof[i] = dmmoab->lidmap[(PetscInt)points[i] - dmmoab->seqstart];
492: PetscFunctionReturn(PETSC_SUCCESS);
493: }
495: /*@C
496: DMMoabGetVertexDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
497: array of locally owned MOAB mesh vertices.
499: Not Collective
501: Input Parameters:
502: . dm - the discretization manager object
504: Output Parameter:
505: . dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
507: Level: intermediate
509: Note:
510: It's utility is when performing Finite-Difference type calculations where vertex traversal is
511: faster than element-wise assembly that is typically done in FEM calculations.
513: .seealso: `DMMoabGetVertexDofsBlockedLocal()`, `DMMoabGetDofsBlocked()`, `DMMoabGetDofsBlockedLocal()`
514: @*/
515: PetscErrorCode DMMoabGetVertexDofsBlocked(DM dm, PetscInt **dof)
516: {
517: DM_Moab *dmmoab;
519: PetscFunctionBegin;
521: dmmoab = (DM_Moab *)dm->data;
523: *dof = dmmoab->gidmap;
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: /*@C
528: DMMoabGetVertexDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
529: array of locally owned MOAB mesh vertices.
531: Not Collective
533: Input Parameters:
534: . dm - the discretization manager object
536: Output Parameter:
537: . dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
539: Level: intermediate
541: Note:
542: It's utility is when performing Finite-Difference type calculations where vertex traversal is
543: faster than element-wise assembly that is typically done in FEM calculations.
545: .seealso: `DMMoabGetVertexDofsBlocked()`, `DMMoabGetDofsBlocked()`, `DMMoabGetDofsBlockedLocal()`
546: @*/
547: PetscErrorCode DMMoabGetVertexDofsBlockedLocal(DM dm, PetscInt **dof)
548: {
549: DM_Moab *dmmoab;
551: PetscFunctionBegin;
553: PetscAssertPointer(dof, 2);
554: dmmoab = (DM_Moab *)dm->data;
556: *dof = dmmoab->lidmap;
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }