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;
30: dmmoab = (DM_Moab *)(dm)->data;
34: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
35: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
36: MBERRNM(merr);
38: DMMoabGetVecTag(fvec, &vtag);
40: merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
41: if (!tag_name.length() && merr != moab::MB_SUCCESS) {
42: VecGetArrayRead(fvec, &varray);
43: /* use the entity handle and the Dof index to set the right value */
44: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)varray);
45: MBERRNM(merr);
46: VecRestoreArrayRead(fvec, &varray);
47: } else {
48: PetscMalloc1(dmmoab->nloc, &farray);
49: /* we are using a MOAB Vec - directly copy the tag data to new one */
50: merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void *)farray);
51: MBERRNM(merr);
52: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)farray);
53: MBERRNM(merr);
54: /* make sure the parallel exchange for ghosts are done appropriately */
55: PetscFree(farray);
56: }
57: #ifdef MOAB_HAVE_MPI
58: merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vowned);
59: MBERRNM(merr);
60: #endif
61: return 0;
62: }
64: /*@C
65: DMMoabSetGlobalFieldVector - Sets the vector reference that represents the global solution associated
66: with all fields (components) managed by DM.
67: A useful utility when updating the DM solution after a solve, to be serialized with the mesh for
68: checkpointing purposes.
70: Not Collective
72: Input Parameters:
73: + dm - the discretization manager object
74: - fvec - the global Vector solution corresponding to all the fields managed by DM
76: Level: intermediate
78: .seealso: `DMMoabGetFieldName()`, `DMMoabSetFieldVector()`
79: @*/
80: PetscErrorCode DMMoabSetGlobalFieldVector(DM dm, Vec fvec)
81: {
82: DM_Moab *dmmoab;
83: moab::Tag vtag, ntag;
84: const PetscScalar *rarray;
85: PetscScalar *varray, *farray;
86: moab::ErrorCode merr;
87: PetscInt i, ifield;
88: std::string tag_name;
89: moab::Range::iterator iter;
92: dmmoab = (DM_Moab *)(dm)->data;
94: /* get the Tag corresponding to the global vector - possible that there is no tag associated.. */
95: DMMoabGetVecTag(fvec, &vtag);
96: merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
97: PetscMalloc1(dmmoab->nloc, &farray);
98: if (!tag_name.length() && merr != moab::MB_SUCCESS) {
99: /* not a MOAB vector - use VecGetSubVector to get the parts as needed */
100: VecGetArrayRead(fvec, &rarray);
101: for (ifield = 0; ifield < dmmoab->numFields; ++ifield) {
102: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
103: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
104: MBERRNM(merr);
106: for (i = 0; i < dmmoab->nloc; i++) farray[i] = (dmmoab->bs == 1 ? rarray[ifield * dmmoab->nloc + i] : rarray[i * dmmoab->numFields + ifield]);
108: /* use the entity handle and the Dof index to set the right value */
109: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)farray);
110: MBERRNM(merr);
111: }
112: VecRestoreArrayRead(fvec, &rarray);
113: } else {
114: PetscMalloc1(dmmoab->nloc * dmmoab->numFields, &varray);
116: /* we are using a MOAB Vec - directly copy the tag data to new one */
117: merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void *)varray);
118: MBERRNM(merr);
119: for (ifield = 0; ifield < dmmoab->numFields; ++ifield) {
120: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
121: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
122: MBERRNM(merr);
124: /* we are using a MOAB Vec - directly copy the tag data to new one */
125: for (i = 0; i < dmmoab->nloc; i++) farray[i] = (dmmoab->bs == 1 ? varray[ifield * dmmoab->nloc + i] : varray[i * dmmoab->numFields + ifield]);
127: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)farray);
128: MBERRNM(merr);
130: #ifdef MOAB_HAVE_MPI
131: /* make sure the parallel exchange for ghosts are done appropriately */
132: merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vlocal);
133: MBERRNM(merr);
134: #endif
135: }
136: PetscFree(varray);
137: }
138: PetscFree(farray);
139: return 0;
140: }
142: /*@C
143: DMMoabSetFieldNames - Sets the number of fields and their names to be managed by the DM
145: Not Collective
147: Input Parameters:
148: + dm - the discretization manager object
149: . numFields - the total number of fields
150: - fields - the array containing the names of each field (component); Can be NULL.
152: Level: intermediate
154: .seealso: `DMMoabGetFieldName()`, `DMMoabSetFieldName()`
155: @*/
156: PetscErrorCode DMMoabSetFieldNames(DM dm, PetscInt numFields, const char *fields[])
157: {
158: PetscInt i;
159: DM_Moab *dmmoab;
162: dmmoab = (DM_Moab *)(dm)->data;
164: /* first deallocate the existing field structure */
165: if (dmmoab->fieldNames) {
166: for (i = 0; i < dmmoab->numFields; i++) PetscFree(dmmoab->fieldNames[i]);
167: PetscFree(dmmoab->fieldNames);
168: }
170: /* now re-allocate and assign field names */
171: dmmoab->numFields = numFields;
172: PetscMalloc1(numFields, &dmmoab->fieldNames);
173: if (fields) {
174: for (i = 0; i < dmmoab->numFields; i++) PetscStrallocpy(fields[i], (char **)&dmmoab->fieldNames[i]);
175: }
176: DMSetNumFields(dm, numFields);
177: return 0;
178: }
180: /*@C
181: DMMoabGetFieldName - Gets the names of individual field components in multicomponent
182: vectors associated with a DMDA.
184: Not Collective
186: Input Parameters:
187: + dm - the discretization manager object
188: - field - field number for the DMMoab (0, 1, ... dof-1), where dof indicates the
189: number of degrees of freedom per node within the DMMoab
191: Output Parameter:
192: . fieldName - the name of the field (component)
194: Level: intermediate
196: .seealso: `DMMoabSetFieldName()`, `DMMoabSetFields()`
197: @*/
198: PetscErrorCode DMMoabGetFieldName(DM dm, PetscInt field, const char **fieldName)
199: {
200: DM_Moab *dmmoab;
203: dmmoab = (DM_Moab *)(dm)->data;
206: *fieldName = dmmoab->fieldNames[field];
207: return 0;
208: }
210: /*@C
211: DMMoabSetFieldName - Sets the name of a field (component) managed by the DM
213: Not Collective
215: Input Parameters:
216: + dm - the discretization manager object
217: . field - the field number
218: - fieldName - the field (component) name
220: Level: intermediate
221: Notes:
222: Can only be called after DMMoabSetFields supplied with correct numFields
224: .seealso: `DMMoabGetFieldName()`, `DMMoabSetFields()`
225: @*/
226: PetscErrorCode DMMoabSetFieldName(DM dm, PetscInt field, const char *fieldName)
227: {
228: DM_Moab *dmmoab;
233: dmmoab = (DM_Moab *)(dm)->data;
236: if (dmmoab->fieldNames[field]) PetscFree(dmmoab->fieldNames[field]);
237: PetscStrallocpy(fieldName, (char **)&dmmoab->fieldNames[field]);
238: return 0;
239: }
241: /*@C
242: DMMoabGetFieldDof - Gets the global degree-of-freedom of a field (component) defined on a
243: particular MOAB EntityHandle.
245: Not Collective
247: Input Parameters:
248: + dm - the discretization manager object
249: . point - the MOAB EntityHandle container which holds the field degree-of-freedom values
250: - field - the field (component) index
252: Output Parameter:
253: . dof - the global degree-of-freedom index corresponding to the field in the discrete representation (Vec, Mat)
255: Level: beginner
257: .seealso: `DMMoabGetFieldDofs()`, `DMMoabGetFieldDofsLocal()`
258: @*/
259: PetscErrorCode DMMoabGetFieldDof(DM dm, moab::EntityHandle point, PetscInt field, PetscInt *dof)
260: {
261: DM_Moab *dmmoab;
264: dmmoab = (DM_Moab *)(dm)->data;
266: *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);
267: return 0;
268: }
270: /*@C
271: DMMoabGetFieldDofs - Gets the global degree-of-freedom of a field (component) defined on an
272: array of MOAB EntityHandles.
274: Not Collective
276: Input Parameters:
277: + dm - the discretization manager object
278: . npoints - the total number of Entities in the points array
279: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
280: - field - the field (component) index
282: Output Parameter:
283: . dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
285: Level: intermediate
287: .seealso: `DMMoabGetFieldDof()`, `DMMoabGetFieldDofsLocal()`
288: @*/
289: PetscErrorCode DMMoabGetFieldDofs(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt field, PetscInt *dof)
290: {
291: PetscInt i;
292: DM_Moab *dmmoab;
296: dmmoab = (DM_Moab *)(dm)->data;
298: if (!dof) PetscMalloc1(npoints, &dof);
300: /* compute the DOF based on local blocking in the fields */
301: /* We also assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
302: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
303: for (i = 0; i < npoints; ++i)
304: 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);
305: return 0;
306: }
308: /*@C
309: DMMoabGetFieldDofsLocal - Gets the local degrees-of-freedom of a field (component) defined on an
310: array of MOAB EntityHandles.
312: Not Collective
314: Input Parameters:
315: + dm - the discretization manager object
316: . npoints - the total number of Entities in the points array
317: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
318: - field - the field (component) index
320: Output Parameter:
321: . dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
323: Level: intermediate
325: .seealso: `DMMoabGetFieldDof()`, `DMMoabGetFieldDofs()`
326: @*/
327: PetscErrorCode DMMoabGetFieldDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt field, PetscInt *dof)
328: {
329: PetscInt i;
330: DM_Moab *dmmoab;
334: dmmoab = (DM_Moab *)(dm)->data;
336: if (!dof) PetscMalloc1(npoints, &dof);
338: /* compute the DOF based on local blocking in the fields */
339: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
340: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
341: for (i = 0; i < npoints; ++i) {
342: 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);
343: }
344: return 0;
345: }
347: /*@C
348: DMMoabGetDofs - Gets the global degree-of-freedom for all fields (components) defined on an
349: array of MOAB EntityHandles.
351: Not Collective
353: Input Parameters:
354: + dm - the discretization manager object
355: . npoints - the total number of Entities in the points array
356: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
358: Output Parameter:
359: . dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
361: Level: intermediate
363: .seealso: `DMMoabGetFieldDofs()`, `DMMoabGetDofsLocal()`, `DMMoabGetDofsBlocked()`
364: @*/
365: PetscErrorCode DMMoabGetDofs(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
366: {
367: PetscInt i, field, offset;
368: DM_Moab *dmmoab;
372: dmmoab = (DM_Moab *)(dm)->data;
374: if (!dof) PetscMalloc1(dmmoab->numFields * npoints, &dof);
376: /* compute the DOF based on local blocking in the fields */
377: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
378: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
379: for (field = 0; field < dmmoab->numFields; ++field) {
380: offset = field * dmmoab->n;
381: for (i = 0; i < npoints; ++i)
382: 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);
383: }
384: return 0;
385: }
387: /*@C
388: DMMoabGetDofsLocal - Gets the local degree-of-freedom for all fields (components) defined on an
389: array of MOAB EntityHandles.
391: Not Collective
393: Input Parameters:
394: + dm - the discretization manager object
395: . npoints - the total number of Entities in the points array
396: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
398: Output Parameter:
399: . dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
401: Level: intermediate
403: .seealso: `DMMoabGetFieldDofs()`, `DMMoabGetDofs()`, `DMMoabGetDofsBlocked()`
404: @*/
405: PetscErrorCode DMMoabGetDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
406: {
407: PetscInt i, field, offset;
408: DM_Moab *dmmoab;
412: dmmoab = (DM_Moab *)(dm)->data;
414: if (!dof) PetscMalloc1(dmmoab->numFields * npoints, &dof);
416: /* compute the DOF based on local blocking in the fields */
417: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
418: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
419: for (field = 0; field < dmmoab->numFields; ++field) {
420: offset = field * dmmoab->n;
421: for (i = 0; i < npoints; ++i)
422: 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);
423: }
424: return 0;
425: }
427: /*@C
428: DMMoabGetDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
429: array of MOAB EntityHandles. It is useful when performing Blocked(Get/Set) methods in computation
430: of element residuals and assembly of the discrete systems when all fields are co-located.
432: Not Collective
434: Input Parameters:
435: + dm - the discretization manager object
436: . npoints - the total number of Entities in the points array
437: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
439: Output Parameter:
440: . dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat)
442: Level: intermediate
444: .seealso: `DMMoabGetDofsLocal()`, `DMMoabGetDofs()`, `DMMoabGetDofsBlockedLocal()`
445: @*/
446: PetscErrorCode DMMoabGetDofsBlocked(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
447: {
448: PetscInt i;
449: DM_Moab *dmmoab;
453: dmmoab = (DM_Moab *)(dm)->data;
455: if (!dof) PetscMalloc1(npoints, &dof);
457: for (i = 0; i < npoints; ++i) dof[i] = dmmoab->gidmap[(PetscInt)points[i] - dmmoab->seqstart];
458: return 0;
459: }
461: /*@C
462: DMMoabGetDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
463: array of MOAB EntityHandles. It is useful when performing local Blocked(Get/Set) methods in computation
464: of element residuals and assembly of the discrete systems when all fields are co-located.
466: Not Collective
468: Input Parameters:
469: + dm - the discretization manager object
470: . npoints - the total number of Entities in the points array
471: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
473: Output Parameter:
474: . dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat)
476: Level: intermediate
478: .seealso: `DMMoabGetDofsLocal()`, `DMMoabGetDofs()`, `DMMoabGetDofsBlockedLocal()`
479: @*/
480: PetscErrorCode DMMoabGetDofsBlockedLocal(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
481: {
482: PetscInt i;
483: DM_Moab *dmmoab;
487: dmmoab = (DM_Moab *)(dm)->data;
489: if (!dof) PetscMalloc1(npoints, &dof);
491: for (i = 0; i < npoints; ++i) dof[i] = dmmoab->lidmap[(PetscInt)points[i] - dmmoab->seqstart];
492: return 0;
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. It's utility is when performing Finite-Difference type calculations
498: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
500: Not Collective
502: Input Parameters:
503: . dm - the discretization manager object
505: Output Parameter:
506: . dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
508: Level: intermediate
510: .seealso: `DMMoabGetVertexDofsBlockedLocal()`, `DMMoabGetDofsBlocked()`, `DMMoabGetDofsBlockedLocal()`
511: @*/
512: PetscErrorCode DMMoabGetVertexDofsBlocked(DM dm, PetscInt **dof)
513: {
514: DM_Moab *dmmoab;
517: dmmoab = (DM_Moab *)(dm)->data;
519: *dof = dmmoab->gidmap;
520: return 0;
521: }
523: /*@C
524: DMMoabGetVertexDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
525: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
526: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
528: Not Collective
530: Input Parameters:
531: . dm - the discretization manager object
533: Output Parameter:
534: . dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
536: Level: intermediate
538: .seealso: `DMMoabGetVertexDofsBlocked()`, `DMMoabGetDofsBlocked()`, `DMMoabGetDofsBlockedLocal()`
539: @*/
540: PetscErrorCode DMMoabGetVertexDofsBlockedLocal(DM dm, PetscInt **dof)
541: {
542: DM_Moab *dmmoab;
546: dmmoab = (DM_Moab *)(dm)->data;
548: *dof = dmmoab->lidmap;
549: return 0;
550: }