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: PetscErrorCode ierr;
28: std::string tag_name;
32: dmmoab = (DM_Moab*)(dm)->data;
34: if ((ifield < 0) || (ifield >= dmmoab->numFields)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The field %d should be positive and less than %d.", ifield, dmmoab->numFields);
36: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
37: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag,
38: moab::MB_TAG_DENSE | moab::MB_TAG_CREAT); MBERRNM(merr);
40: DMMoabGetVecTag(fvec, &vtag);
42: merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
43: if (!tag_name.length() && merr != moab::MB_SUCCESS) {
44: VecGetArrayRead(fvec, &varray);
45: /* use the entity handle and the Dof index to set the right value */
46: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)varray); MBERRNM(merr);
47: VecRestoreArrayRead(fvec, &varray);
48: }
49: else {
50: PetscMalloc1(dmmoab->nloc, &farray);
51: /* we are using a MOAB Vec - directly copy the tag data to new one */
52: merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void*)farray); MBERRNM(merr);
53: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray); 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); MBERRNM(merr);
59: #endif
60: return(0);
61: }
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: PetscErrorCode ierr;
88: PetscInt i, ifield;
89: std::string tag_name;
90: moab::Range::iterator iter;
94: dmmoab = (DM_Moab*)(dm)->data;
96: /* get the Tag corresponding to the global vector - possible that there is no tag associated.. */
97: DMMoabGetVecTag(fvec, &vtag);
98: merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
99: 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: VecGetArrayRead(fvec, &rarray);
103: for (ifield = 0; ifield < dmmoab->numFields; ++ifield) {
105: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
106: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag,
107: moab::MB_TAG_DENSE | moab::MB_TAG_CREAT); MBERRNM(merr);
109: for (i = 0; i < dmmoab->nloc; i++) {
110: farray[i] = (dmmoab->bs == 1 ? rarray[ifield * dmmoab->nloc + i] : rarray[i * dmmoab->numFields + ifield]);
111: }
113: /* use the entity handle and the Dof index to set the right value */
114: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray); MBERRNM(merr);
115: }
116: VecRestoreArrayRead(fvec, &rarray);
117: }
118: else {
119: PetscMalloc1(dmmoab->nloc * dmmoab->numFields, &varray);
121: /* we are using a MOAB Vec - directly copy the tag data to new one */
122: merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void*)varray); MBERRNM(merr);
123: for (ifield = 0; ifield < dmmoab->numFields; ++ifield) {
125: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
126: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag,
127: moab::MB_TAG_DENSE | moab::MB_TAG_CREAT); MBERRNM(merr);
129: /* we are using a MOAB Vec - directly copy the tag data to new one */
130: for (i = 0; i < dmmoab->nloc; i++) {
131: farray[i] = (dmmoab->bs == 1 ? varray[ifield * dmmoab->nloc + i] : varray[i * dmmoab->numFields + ifield]);
132: }
134: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray); MBERRNM(merr);
136: #ifdef MOAB_HAVE_MPI
137: /* make sure the parallel exchange for ghosts are done appropriately */
138: merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vlocal); MBERRNM(merr);
139: #endif
140: }
141: PetscFree(varray);
142: }
143: PetscFree(farray);
144: return(0);
145: }
148: /*@C
149: DMMoabSetFieldNames - Sets the number of fields and their names to be managed by the DM
151: Not Collective
153: Input Parameters:
154: + dm - the discretization manager object
155: . numFields - the total number of fields
156: - fields - the array containing the names of each field (component); Can be NULL.
158: Level: intermediate
160: .seealso: DMMoabGetFieldName(), DMMoabSetFieldName()
161: @*/
162: PetscErrorCode DMMoabSetFieldNames(DM dm, PetscInt numFields, const char* fields[])
163: {
165: PetscInt i;
166: DM_Moab *dmmoab;
170: dmmoab = (DM_Moab*)(dm)->data;
172: /* first deallocate the existing field structure */
173: if (dmmoab->fieldNames) {
174: for (i = 0; i < dmmoab->numFields; i++) {
175: PetscFree(dmmoab->fieldNames[i]);
176: }
177: PetscFree(dmmoab->fieldNames);
178: }
180: /* now re-allocate and assign field names */
181: dmmoab->numFields = numFields;
182: PetscMalloc1(numFields, &dmmoab->fieldNames);
183: if (fields) {
184: for (i = 0; i < dmmoab->numFields; i++) {
185: PetscStrallocpy(fields[i], (char**) &dmmoab->fieldNames[i]);
186: }
187: }
188: DMSetNumFields(dm, numFields);
189: return(0);
190: }
193: /*@C
194: DMMoabGetFieldName - Gets the names of individual field components in multicomponent
195: vectors associated with a DMDA.
197: Not Collective
199: Input Parameter:
200: + dm - the discretization manager object
201: - field - field number for the DMMoab (0, 1, ... dof-1), where dof indicates the
202: number of degrees of freedom per node within the DMMoab
204: Output Parameter:
205: . fieldName - the name of the field (component)
207: Level: intermediate
209: .seealso: DMMoabSetFieldName(), DMMoabSetFields()
210: @*/
211: PetscErrorCode DMMoabGetFieldName(DM dm, PetscInt field, const char **fieldName)
212: {
213: DM_Moab *dmmoab;
217: dmmoab = (DM_Moab*)(dm)->data;
218: if ((field < 0) || (field >= dmmoab->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "DM field %d should be in [%d, %d)", field, 0, dmmoab->numFields);
220: *fieldName = dmmoab->fieldNames[field];
221: return(0);
222: }
225: /*@C
226: DMMoabSetFieldName - Sets the name of a field (component) managed by the DM
228: Not Collective
230: Input Parameters:
231: + dm - the discretization manager object
232: . field - the field number
233: - fieldName - the field (component) name
235: Level: intermediate
236: Notes:
237: Can only be called after DMMoabSetFields supplied with correct numFields
239: .seealso: DMMoabGetFieldName(), DMMoabSetFields()
240: @*/
241: PetscErrorCode DMMoabSetFieldName(DM dm, PetscInt field, const char *fieldName)
242: {
244: DM_Moab *dmmoab;
250: dmmoab = (DM_Moab*)(dm)->data;
251: if ((field < 0) || (field >= dmmoab->numFields)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "DM field %d should be in [%d, %d)", field, 0, dmmoab->numFields);
253: if (dmmoab->fieldNames[field]) {
254: PetscFree(dmmoab->fieldNames[field]);
255: }
256: PetscStrallocpy(fieldName, (char**) &dmmoab->fieldNames[field]);
257: return(0);
258: }
261: /*@C
262: DMMoabGetFieldDof - Gets the global degree-of-freedom of a field (component) defined on a
263: particular MOAB EntityHandle.
265: Not Collective
267: Input Parameters:
268: + dm - the discretization manager object
269: . point - the MOAB EntityHandle container which holds the field degree-of-freedom values
270: - field - the field (component) index
272: Output Parameter:
273: . dof - the global degree-of-freedom index corresponding to the field in the discrete representation (Vec, Mat)
275: Level: beginner
277: .seealso: DMMoabGetFieldDofs(), DMMoabGetFieldDofsLocal()
278: @*/
279: PetscErrorCode DMMoabGetFieldDof(DM dm, moab::EntityHandle point, PetscInt field, PetscInt* dof)
280: {
281: DM_Moab *dmmoab;
285: dmmoab = (DM_Moab*)(dm)->data;
287: *dof = (dmmoab->bs == 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(point) - dmmoab->seqstart] + field * dmmoab->n :
288: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(point) - dmmoab->seqstart] * dmmoab->numFields + field);
289: return(0);
290: }
293: /*@C
294: DMMoabGetFieldDofs - Gets the global degree-of-freedom of a field (component) defined on an
295: array of MOAB EntityHandles.
297: Not Collective
299: Input Parameters:
300: + dm - the discretization manager object
301: . npoints - the total number of Entities in the points array
302: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
303: - field - the field (component) index
305: Output Parameter:
306: . dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
308: Level: intermediate
310: .seealso: DMMoabGetFieldDof(), DMMoabGetFieldDofsLocal()
311: @*/
312: PetscErrorCode DMMoabGetFieldDofs(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt field, PetscInt* dof)
313: {
314: PetscInt i;
315: PetscErrorCode ierr;
316: DM_Moab *dmmoab;
321: dmmoab = (DM_Moab*)(dm)->data;
323: if (!dof) {
324: PetscMalloc1(npoints, &dof);
325: }
327: /* compute the DOF based on local blocking in the fields */
328: /* We also assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
329: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
330: for (i = 0; i < npoints; ++i)
331: dof[i] = (dmmoab->bs == 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + field * dmmoab->n :
332: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field);
333: return(0);
334: }
337: /*@C
338: DMMoabGetFieldDofsLocal - Gets the local degrees-of-freedom of a field (component) defined on an
339: array of MOAB EntityHandles.
341: Not Collective
343: Input Parameters:
344: + dm - the discretization manager object
345: . npoints - the total number of Entities in the points array
346: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
347: - field - the field (component) index
349: Output Parameter:
350: . dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
352: Level: intermediate
354: .seealso: DMMoabGetFieldDof(), DMMoabGetFieldDofs()
355: @*/
356: PetscErrorCode DMMoabGetFieldDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt field, PetscInt* dof)
357: {
358: PetscInt i;
359: PetscErrorCode ierr;
360: DM_Moab *dmmoab;
365: dmmoab = (DM_Moab*)(dm)->data;
367: if (!dof) {
368: PetscMalloc1(npoints, &dof);
369: }
371: /* compute the DOF based on local blocking in the fields */
372: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
373: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
374: for (i = 0; i < npoints; ++i) {
375: dof[i] = (dmmoab->bs > 1 ? dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
376: dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + field * dmmoab->n);
377: }
378: return(0);
379: }
382: /*@C
383: DMMoabGetDofs - Gets the global degree-of-freedom for all fields (components) defined on an
384: array of MOAB EntityHandles.
386: Not Collective
388: Input Parameters:
389: + dm - the discretization manager object
390: . npoints - the total number of Entities in the points array
391: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
393: Output Parameter:
394: . dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
396: Level: intermediate
398: .seealso: DMMoabGetFieldDofs(), DMMoabGetDofsLocal(), DMMoabGetDofsBlocked()
399: @*/
400: PetscErrorCode DMMoabGetDofs(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
401: {
402: PetscInt i, field, offset;
403: PetscErrorCode ierr;
404: DM_Moab *dmmoab;
409: dmmoab = (DM_Moab*)(dm)->data;
411: if (!dof) {
412: PetscMalloc1(dmmoab->numFields * npoints, &dof);
413: }
415: /* compute the DOF based on local blocking in the fields */
416: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
417: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
418: for (field = 0; field < dmmoab->numFields; ++field) {
419: offset = field * dmmoab->n;
420: for (i = 0; i < npoints; ++i)
421: dof[i * dmmoab->numFields + field] = (dmmoab->bs > 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
422: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + offset);
423: }
424: return(0);
425: }
427: /*@C
428: DMMoabGetDofsLocal - Gets the local degree-of-freedom for all fields (components) defined on an
429: array of MOAB EntityHandles.
431: Not Collective
433: Input Parameters:
434: + dm - the discretization manager object
435: . npoints - the total number of Entities in the points array
436: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
438: Output Parameter:
439: . dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
441: Level: intermediate
443: .seealso: DMMoabGetFieldDofs(), DMMoabGetDofs(), DMMoabGetDofsBlocked()
444: @*/
445: PetscErrorCode DMMoabGetDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
446: {
447: PetscInt i, field, offset;
448: PetscErrorCode ierr;
449: DM_Moab *dmmoab;
454: dmmoab = (DM_Moab*)(dm)->data;
456: if (!dof) {
457: PetscMalloc1(dmmoab->numFields * npoints, &dof);
458: }
460: /* compute the DOF based on local blocking in the fields */
461: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
462: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
463: for (field = 0; field < dmmoab->numFields; ++field) {
464: offset = field * dmmoab->n;
465: for (i = 0; i < npoints; ++i)
466: dof[i * dmmoab->numFields + field] = (dmmoab->bs > 1 ? dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
467: dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + offset);
468: }
469: return(0);
470: }
473: /*@C
474: DMMoabGetDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
475: array of MOAB EntityHandles. It is useful when performing Blocked(Get/Set) methods in computation
476: of element residuals and assembly of the discrete systems when all fields are co-located.
478: Not Collective
480: Input Parameters:
481: + dm - the discretization manager object
482: . npoints - the total number of Entities in the points array
483: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
485: Output Parameter:
486: . dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat)
488: Level: intermediate
490: .seealso: DMMoabGetDofsLocal(), DMMoabGetDofs(), DMMoabGetDofsBlockedLocal()
491: @*/
492: PetscErrorCode DMMoabGetDofsBlocked(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
493: {
494: PetscInt i;
495: DM_Moab *dmmoab;
496: PetscErrorCode ierr;
501: dmmoab = (DM_Moab*)(dm)->data;
503: if (!dof) {
504: PetscMalloc1(npoints, &dof);
505: }
507: for (i = 0; i < npoints; ++i) {
508: dof[i] = dmmoab->gidmap[(PetscInt)points[i] - dmmoab->seqstart];
509: }
510: return(0);
511: }
514: /*@C
515: DMMoabGetDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
516: array of MOAB EntityHandles. It is useful when performing local Blocked(Get/Set) methods in computation
517: of element residuals and assembly of the discrete systems when all fields are co-located.
519: Not Collective
521: Input Parameters:
522: + dm - the discretization manager object
523: . npoints - the total number of Entities in the points array
524: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
526: Output Parameter:
527: . dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat)
529: Level: intermediate
531: .seealso: DMMoabGetDofsLocal(), DMMoabGetDofs(), DMMoabGetDofsBlockedLocal()
532: @*/
533: PetscErrorCode DMMoabGetDofsBlockedLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
534: {
535: PetscInt i;
536: DM_Moab *dmmoab;
537: PetscErrorCode ierr;
542: dmmoab = (DM_Moab*)(dm)->data;
544: if (!dof) {
545: PetscMalloc1(npoints, &dof);
546: }
548: for (i = 0; i < npoints; ++i)
549: dof[i] = dmmoab->lidmap[(PetscInt)points[i] - dmmoab->seqstart];
550: return(0);
551: }
554: /*@C
555: DMMoabGetVertexDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
556: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
557: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
559: Not Collective
561: Input Parameters:
562: . dm - the discretization manager object
564: Output Parameter:
565: . dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
567: Level: intermediate
569: .seealso: DMMoabGetVertexDofsBlockedLocal(), DMMoabGetDofsBlocked(), DMMoabGetDofsBlockedLocal()
570: @*/
571: PetscErrorCode DMMoabGetVertexDofsBlocked(DM dm, PetscInt** dof)
572: {
573: DM_Moab *dmmoab;
577: dmmoab = (DM_Moab*)(dm)->data;
579: *dof = dmmoab->gidmap;
580: return(0);
581: }
584: /*@C
585: DMMoabGetVertexDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
586: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
587: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
589: Not Collective
591: Input Parameters:
592: . dm - the discretization manager object
594: Output Parameter:
595: . dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
597: Level: intermediate
599: .seealso: DMMoabGetVertexDofsBlocked(), DMMoabGetDofsBlocked(), DMMoabGetDofsBlockedLocal()
600: @*/
601: PetscErrorCode DMMoabGetVertexDofsBlockedLocal(DM dm, PetscInt** dof)
602: {
603: DM_Moab *dmmoab;
608: dmmoab = (DM_Moab*)(dm)->data;
610: *dof = dmmoab->lidmap;
611: return(0);
612: }