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,
36: moab::MB_TAG_DENSE | moab::MB_TAG_CREAT); 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); MBERRNM(merr);
45: VecRestoreArrayRead(fvec, &varray);
46: }
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); MBERRNM(merr);
51: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray); MBERRNM(merr);
52: /* make sure the parallel exchange for ghosts are done appropriately */
53: PetscFree(farray);
54: }
55: #ifdef MOAB_HAVE_MPI
56: merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vowned); MBERRNM(merr);
57: #endif
58: return 0;
59: }
61: /*@C
62: DMMoabSetGlobalFieldVector - Sets the vector reference that represents the global solution associated
63: with all fields (components) managed by DM.
64: A useful utility when updating the DM solution after a solve, to be serialized with the mesh for
65: checkpointing purposes.
67: Not Collective
69: Input Parameters:
70: + dm - the discretization manager object
71: - fvec - the global Vector solution corresponding to all the fields managed by DM
73: Level: intermediate
75: .seealso: DMMoabGetFieldName(), DMMoabSetFieldVector()
76: @*/
77: PetscErrorCode DMMoabSetGlobalFieldVector(DM dm, Vec fvec)
78: {
79: DM_Moab *dmmoab;
80: moab::Tag vtag, ntag;
81: const PetscScalar *rarray;
82: PetscScalar *varray, *farray;
83: moab::ErrorCode merr;
84: PetscInt i, ifield;
85: std::string tag_name;
86: moab::Range::iterator iter;
89: dmmoab = (DM_Moab*)(dm)->data;
91: /* get the Tag corresponding to the global vector - possible that there is no tag associated.. */
92: DMMoabGetVecTag(fvec, &vtag);
93: merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
94: PetscMalloc1(dmmoab->nloc, &farray);
95: if (!tag_name.length() && merr != moab::MB_SUCCESS) {
96: /* not a MOAB vector - use VecGetSubVector to get the parts as needed */
97: VecGetArrayRead(fvec, &rarray);
98: for (ifield = 0; ifield < dmmoab->numFields; ++ifield) {
100: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
101: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag,
102: moab::MB_TAG_DENSE | moab::MB_TAG_CREAT); MBERRNM(merr);
104: for (i = 0; i < dmmoab->nloc; i++) {
105: farray[i] = (dmmoab->bs == 1 ? rarray[ifield * dmmoab->nloc + i] : rarray[i * dmmoab->numFields + ifield]);
106: }
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); MBERRNM(merr);
110: }
111: VecRestoreArrayRead(fvec, &rarray);
112: }
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); MBERRNM(merr);
118: 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,
122: moab::MB_TAG_DENSE | moab::MB_TAG_CREAT); 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++) {
126: farray[i] = (dmmoab->bs == 1 ? varray[ifield * dmmoab->nloc + i] : varray[i * dmmoab->numFields + ifield]);
127: }
129: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray); MBERRNM(merr);
131: #ifdef MOAB_HAVE_MPI
132: /* make sure the parallel exchange for ghosts are done appropriately */
133: merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vlocal); 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++) {
167: PetscFree(dmmoab->fieldNames[i]);
168: }
169: PetscFree(dmmoab->fieldNames);
170: }
172: /* now re-allocate and assign field names */
173: dmmoab->numFields = numFields;
174: PetscMalloc1(numFields, &dmmoab->fieldNames);
175: if (fields) {
176: for (i = 0; i < dmmoab->numFields; i++) {
177: PetscStrallocpy(fields[i], (char**) &dmmoab->fieldNames[i]);
178: }
179: }
180: DMSetNumFields(dm, numFields);
181: return 0;
182: }
184: /*@C
185: DMMoabGetFieldName - Gets the names of individual field components in multicomponent
186: vectors associated with a DMDA.
188: Not Collective
190: Input Parameters:
191: + dm - the discretization manager object
192: - field - field number for the DMMoab (0, 1, ... dof-1), where dof indicates the
193: number of degrees of freedom per node within the DMMoab
195: Output Parameter:
196: . fieldName - the name of the field (component)
198: Level: intermediate
200: .seealso: DMMoabSetFieldName(), DMMoabSetFields()
201: @*/
202: PetscErrorCode DMMoabGetFieldName(DM dm, PetscInt field, const char **fieldName)
203: {
204: DM_Moab *dmmoab;
207: dmmoab = (DM_Moab*)(dm)->data;
210: *fieldName = dmmoab->fieldNames[field];
211: return 0;
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;
237: dmmoab = (DM_Moab*)(dm)->data;
240: if (dmmoab->fieldNames[field]) {
241: PetscFree(dmmoab->fieldNames[field]);
242: }
243: PetscStrallocpy(fieldName, (char**) &dmmoab->fieldNames[field]);
244: return 0;
245: }
247: /*@C
248: DMMoabGetFieldDof - Gets the global degree-of-freedom of a field (component) defined on a
249: particular MOAB EntityHandle.
251: Not Collective
253: Input Parameters:
254: + dm - the discretization manager object
255: . point - the MOAB EntityHandle container which holds the field degree-of-freedom values
256: - field - the field (component) index
258: Output Parameter:
259: . dof - the global degree-of-freedom index corresponding to the field in the discrete representation (Vec, Mat)
261: Level: beginner
263: .seealso: DMMoabGetFieldDofs(), DMMoabGetFieldDofsLocal()
264: @*/
265: PetscErrorCode DMMoabGetFieldDof(DM dm, moab::EntityHandle point, PetscInt field, PetscInt* dof)
266: {
267: DM_Moab *dmmoab;
270: dmmoab = (DM_Moab*)(dm)->data;
272: *dof = (dmmoab->bs == 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(point) - dmmoab->seqstart] + field * dmmoab->n :
273: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(point) - dmmoab->seqstart] * dmmoab->numFields + field);
274: return 0;
275: }
277: /*@C
278: DMMoabGetFieldDofs - Gets the global degree-of-freedom of a field (component) defined on an
279: array of MOAB EntityHandles.
281: Not Collective
283: Input Parameters:
284: + dm - the discretization manager object
285: . npoints - the total number of Entities in the points array
286: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
287: - field - the field (component) index
289: Output Parameter:
290: . dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
292: Level: intermediate
294: .seealso: DMMoabGetFieldDof(), DMMoabGetFieldDofsLocal()
295: @*/
296: PetscErrorCode DMMoabGetFieldDofs(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt field, PetscInt* dof)
297: {
298: PetscInt i;
299: DM_Moab *dmmoab;
303: dmmoab = (DM_Moab*)(dm)->data;
305: if (!dof) {
306: PetscMalloc1(npoints, &dof);
307: }
309: /* compute the DOF based on local blocking in the fields */
310: /* We also assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
311: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
312: for (i = 0; i < npoints; ++i)
313: dof[i] = (dmmoab->bs == 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + field * dmmoab->n :
314: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field);
315: return 0;
316: }
318: /*@C
319: DMMoabGetFieldDofsLocal - Gets the local degrees-of-freedom of a field (component) defined on an
320: array of MOAB EntityHandles.
322: Not Collective
324: Input Parameters:
325: + dm - the discretization manager object
326: . npoints - the total number of Entities in the points array
327: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
328: - field - the field (component) index
330: Output Parameter:
331: . dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
333: Level: intermediate
335: .seealso: DMMoabGetFieldDof(), DMMoabGetFieldDofs()
336: @*/
337: PetscErrorCode DMMoabGetFieldDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt field, PetscInt* dof)
338: {
339: PetscInt i;
340: DM_Moab *dmmoab;
344: dmmoab = (DM_Moab*)(dm)->data;
346: if (!dof) {
347: PetscMalloc1(npoints, &dof);
348: }
350: /* compute the DOF based on local blocking in the fields */
351: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
352: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
353: for (i = 0; i < npoints; ++i) {
354: dof[i] = (dmmoab->bs > 1 ? dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
355: dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + field * dmmoab->n);
356: }
357: return 0;
358: }
360: /*@C
361: DMMoabGetDofs - Gets the global degree-of-freedom for all fields (components) defined on an
362: array of MOAB EntityHandles.
364: Not Collective
366: Input Parameters:
367: + dm - the discretization manager object
368: . npoints - the total number of Entities in the points array
369: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
371: Output Parameter:
372: . dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
374: Level: intermediate
376: .seealso: DMMoabGetFieldDofs(), DMMoabGetDofsLocal(), DMMoabGetDofsBlocked()
377: @*/
378: PetscErrorCode DMMoabGetDofs(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
379: {
380: PetscInt i, field, offset;
381: DM_Moab *dmmoab;
385: dmmoab = (DM_Moab*)(dm)->data;
387: if (!dof) {
388: PetscMalloc1(dmmoab->numFields * npoints, &dof);
389: }
391: /* compute the DOF based on local blocking in the fields */
392: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
393: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
394: for (field = 0; field < dmmoab->numFields; ++field) {
395: offset = field * dmmoab->n;
396: for (i = 0; i < npoints; ++i)
397: dof[i * dmmoab->numFields + field] = (dmmoab->bs > 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
398: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + offset);
399: }
400: return 0;
401: }
403: /*@C
404: DMMoabGetDofsLocal - Gets the local degree-of-freedom for all fields (components) defined on an
405: array of MOAB EntityHandles.
407: Not Collective
409: Input Parameters:
410: + dm - the discretization manager object
411: . npoints - the total number of Entities in the points array
412: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
414: Output Parameter:
415: . dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
417: Level: intermediate
419: .seealso: DMMoabGetFieldDofs(), DMMoabGetDofs(), DMMoabGetDofsBlocked()
420: @*/
421: PetscErrorCode DMMoabGetDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
422: {
423: PetscInt i, field, offset;
424: DM_Moab *dmmoab;
428: dmmoab = (DM_Moab*)(dm)->data;
430: if (!dof) {
431: PetscMalloc1(dmmoab->numFields * npoints, &dof);
432: }
434: /* compute the DOF based on local blocking in the fields */
435: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
436: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
437: for (field = 0; field < dmmoab->numFields; ++field) {
438: offset = field * dmmoab->n;
439: for (i = 0; i < npoints; ++i)
440: dof[i * dmmoab->numFields + field] = (dmmoab->bs > 1 ? dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
441: dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + offset);
442: }
443: return 0;
444: }
446: /*@C
447: DMMoabGetDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
448: array of MOAB EntityHandles. It is useful when performing Blocked(Get/Set) methods in computation
449: of element residuals and assembly of the discrete systems when all fields are co-located.
451: Not Collective
453: Input Parameters:
454: + dm - the discretization manager object
455: . npoints - the total number of Entities in the points array
456: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
458: Output Parameter:
459: . dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat)
461: Level: intermediate
463: .seealso: DMMoabGetDofsLocal(), DMMoabGetDofs(), DMMoabGetDofsBlockedLocal()
464: @*/
465: PetscErrorCode DMMoabGetDofsBlocked(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
466: {
467: PetscInt i;
468: DM_Moab *dmmoab;
472: dmmoab = (DM_Moab*)(dm)->data;
474: if (!dof) {
475: PetscMalloc1(npoints, &dof);
476: }
478: for (i = 0; i < npoints; ++i) {
479: dof[i] = dmmoab->gidmap[(PetscInt)points[i] - dmmoab->seqstart];
480: }
481: return 0;
482: }
484: /*@C
485: DMMoabGetDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
486: array of MOAB EntityHandles. It is useful when performing local Blocked(Get/Set) methods in computation
487: of element residuals and assembly of the discrete systems when all fields are co-located.
489: Not Collective
491: Input Parameters:
492: + dm - the discretization manager object
493: . npoints - the total number of Entities in the points array
494: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
496: Output Parameter:
497: . dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat)
499: Level: intermediate
501: .seealso: DMMoabGetDofsLocal(), DMMoabGetDofs(), DMMoabGetDofsBlockedLocal()
502: @*/
503: PetscErrorCode DMMoabGetDofsBlockedLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
504: {
505: PetscInt i;
506: DM_Moab *dmmoab;
510: dmmoab = (DM_Moab*)(dm)->data;
512: if (!dof) PetscMalloc1(npoints, &dof);
514: for (i = 0; i < npoints; ++i) dof[i] = dmmoab->lidmap[(PetscInt)points[i] - dmmoab->seqstart];
515: return 0;
516: }
518: /*@C
519: DMMoabGetVertexDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
520: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
521: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
523: Not Collective
525: Input Parameters:
526: . dm - the discretization manager object
528: Output Parameter:
529: . dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
531: Level: intermediate
533: .seealso: DMMoabGetVertexDofsBlockedLocal(), DMMoabGetDofsBlocked(), DMMoabGetDofsBlockedLocal()
534: @*/
535: PetscErrorCode DMMoabGetVertexDofsBlocked(DM dm, PetscInt** dof)
536: {
537: DM_Moab *dmmoab;
540: dmmoab = (DM_Moab*)(dm)->data;
542: *dof = dmmoab->gidmap;
543: return 0;
544: }
546: /*@C
547: DMMoabGetVertexDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
548: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
549: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
551: Not Collective
553: Input Parameters:
554: . dm - the discretization manager object
556: Output Parameter:
557: . dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
559: Level: intermediate
561: .seealso: DMMoabGetVertexDofsBlocked(), DMMoabGetDofsBlocked(), DMMoabGetDofsBlockedLocal()
562: @*/
563: PetscErrorCode DMMoabGetVertexDofsBlockedLocal(DM dm, PetscInt** dof)
564: {
565: DM_Moab *dmmoab;
569: dmmoab = (DM_Moab*)(dm)->data;
571: *dof = dmmoab->lidmap;
572: return 0;
573: }