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: }
63: /*@C
64: DMMoabSetGlobalFieldVector - Sets the vector reference that represents the global solution associated
65: with all fields (components) managed by DM.
66: A useful utility when updating the DM solution after a solve, to be serialized with the mesh for
67: checkpointing purposes.
69: Not Collective
71: Input Parameters:
72: + dm - the discretization manager object
73: - fvec - the global Vector solution corresponding to all the fields managed by DM
75: Level: intermediate
77: .seealso: DMMoabGetFieldName(), DMMoabSetFieldVector()
78: @*/
79: PetscErrorCode DMMoabSetGlobalFieldVector(DM dm, Vec fvec)
80: {
81: DM_Moab *dmmoab;
82: moab::Tag vtag, ntag;
83: const PetscScalar *rarray;
84: PetscScalar *varray, *farray;
85: moab::ErrorCode merr;
86: PetscErrorCode ierr;
87: PetscInt i, ifield;
88: std::string tag_name;
89: moab::Range::iterator iter;
93: dmmoab = (DM_Moab*)(dm)->data;
95: /* get the Tag corresponding to the global vector - possible that there is no tag associated.. */
96: DMMoabGetVecTag(fvec, &vtag);
97: merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
98: PetscMalloc1(dmmoab->nloc, &farray);
99: if (!tag_name.length() && merr != moab::MB_SUCCESS) {
100: /* not a MOAB vector - use VecGetSubVector to get the parts as needed */
101: VecGetArrayRead(fvec, &rarray);
102: 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,
106: moab::MB_TAG_DENSE | moab::MB_TAG_CREAT); MBERRNM(merr);
108: for (i = 0; i < dmmoab->nloc; i++) {
109: farray[i] = (dmmoab->bs == 1 ? rarray[ifield * dmmoab->nloc + i] : rarray[i * dmmoab->numFields + ifield]);
110: }
112: /* use the entity handle and the Dof index to set the right value */
113: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray); MBERRNM(merr);
114: }
115: VecRestoreArrayRead(fvec, &rarray);
116: }
117: else {
118: PetscMalloc1(dmmoab->nloc * dmmoab->numFields, &varray);
120: /* we are using a MOAB Vec - directly copy the tag data to new one */
121: merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void*)varray); MBERRNM(merr);
122: for (ifield = 0; ifield < dmmoab->numFields; ++ifield) {
124: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
125: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag,
126: moab::MB_TAG_DENSE | moab::MB_TAG_CREAT); MBERRNM(merr);
128: /* we are using a MOAB Vec - directly copy the tag data to new one */
129: for (i = 0; i < dmmoab->nloc; i++) {
130: farray[i] = (dmmoab->bs == 1 ? varray[ifield * dmmoab->nloc + i] : varray[i * dmmoab->numFields + ifield]);
131: }
133: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray); MBERRNM(merr);
135: #ifdef MOAB_HAVE_MPI
136: /* make sure the parallel exchange for ghosts are done appropriately */
137: merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vlocal); MBERRNM(merr);
138: #endif
139: }
140: PetscFree(varray);
141: }
142: PetscFree(farray);
143: return(0);
144: }
146: /*@C
147: DMMoabSetFieldNames - Sets the number of fields and their names to be managed by the DM
149: Not Collective
151: Input Parameters:
152: + dm - the discretization manager object
153: . numFields - the total number of fields
154: - fields - the array containing the names of each field (component); Can be NULL.
156: Level: intermediate
158: .seealso: DMMoabGetFieldName(), DMMoabSetFieldName()
159: @*/
160: PetscErrorCode DMMoabSetFieldNames(DM dm, PetscInt numFields, const char* fields[])
161: {
163: PetscInt i;
164: DM_Moab *dmmoab;
168: dmmoab = (DM_Moab*)(dm)->data;
170: /* first deallocate the existing field structure */
171: if (dmmoab->fieldNames) {
172: for (i = 0; i < dmmoab->numFields; i++) {
173: PetscFree(dmmoab->fieldNames[i]);
174: }
175: PetscFree(dmmoab->fieldNames);
176: }
178: /* now re-allocate and assign field names */
179: dmmoab->numFields = numFields;
180: PetscMalloc1(numFields, &dmmoab->fieldNames);
181: if (fields) {
182: for (i = 0; i < dmmoab->numFields; i++) {
183: PetscStrallocpy(fields[i], (char**) &dmmoab->fieldNames[i]);
184: }
185: }
186: DMSetNumFields(dm, numFields);
187: return(0);
188: }
190: /*@C
191: DMMoabGetFieldName - Gets the names of individual field components in multicomponent
192: vectors associated with a DMDA.
194: Not Collective
196: Input Parameters:
197: + dm - the discretization manager object
198: - field - field number for the DMMoab (0, 1, ... dof-1), where dof indicates the
199: number of degrees of freedom per node within the DMMoab
201: Output Parameter:
202: . fieldName - the name of the field (component)
204: Level: intermediate
206: .seealso: DMMoabSetFieldName(), DMMoabSetFields()
207: @*/
208: PetscErrorCode DMMoabGetFieldName(DM dm, PetscInt field, const char **fieldName)
209: {
210: DM_Moab *dmmoab;
214: dmmoab = (DM_Moab*)(dm)->data;
215: 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);
217: *fieldName = dmmoab->fieldNames[field];
218: return(0);
219: }
221: /*@C
222: DMMoabSetFieldName - Sets the name of a field (component) managed by the DM
224: Not Collective
226: Input Parameters:
227: + dm - the discretization manager object
228: . field - the field number
229: - fieldName - the field (component) name
231: Level: intermediate
232: Notes:
233: Can only be called after DMMoabSetFields supplied with correct numFields
235: .seealso: DMMoabGetFieldName(), DMMoabSetFields()
236: @*/
237: PetscErrorCode DMMoabSetFieldName(DM dm, PetscInt field, const char *fieldName)
238: {
240: DM_Moab *dmmoab;
246: dmmoab = (DM_Moab*)(dm)->data;
247: 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);
249: if (dmmoab->fieldNames[field]) {
250: PetscFree(dmmoab->fieldNames[field]);
251: }
252: PetscStrallocpy(fieldName, (char**) &dmmoab->fieldNames[field]);
253: return(0);
254: }
256: /*@C
257: DMMoabGetFieldDof - Gets the global degree-of-freedom of a field (component) defined on a
258: particular MOAB EntityHandle.
260: Not Collective
262: Input Parameters:
263: + dm - the discretization manager object
264: . point - the MOAB EntityHandle container which holds the field degree-of-freedom values
265: - field - the field (component) index
267: Output Parameter:
268: . dof - the global degree-of-freedom index corresponding to the field in the discrete representation (Vec, Mat)
270: Level: beginner
272: .seealso: DMMoabGetFieldDofs(), DMMoabGetFieldDofsLocal()
273: @*/
274: PetscErrorCode DMMoabGetFieldDof(DM dm, moab::EntityHandle point, PetscInt field, PetscInt* dof)
275: {
276: DM_Moab *dmmoab;
280: dmmoab = (DM_Moab*)(dm)->data;
282: *dof = (dmmoab->bs == 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(point) - dmmoab->seqstart] + field * dmmoab->n :
283: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(point) - dmmoab->seqstart] * dmmoab->numFields + field);
284: return(0);
285: }
287: /*@C
288: DMMoabGetFieldDofs - Gets the global degree-of-freedom of a field (component) defined on an
289: array of MOAB EntityHandles.
291: Not Collective
293: Input Parameters:
294: + dm - the discretization manager object
295: . npoints - the total number of Entities in the points array
296: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
297: - field - the field (component) index
299: Output Parameter:
300: . dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
302: Level: intermediate
304: .seealso: DMMoabGetFieldDof(), DMMoabGetFieldDofsLocal()
305: @*/
306: PetscErrorCode DMMoabGetFieldDofs(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt field, PetscInt* dof)
307: {
308: PetscInt i;
309: PetscErrorCode ierr;
310: DM_Moab *dmmoab;
315: dmmoab = (DM_Moab*)(dm)->data;
317: if (!dof) {
318: PetscMalloc1(npoints, &dof);
319: }
321: /* compute the DOF based on local blocking in the fields */
322: /* We also assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
323: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
324: for (i = 0; i < npoints; ++i)
325: dof[i] = (dmmoab->bs == 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + field * dmmoab->n :
326: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field);
327: return(0);
328: }
330: /*@C
331: DMMoabGetFieldDofsLocal - Gets the local degrees-of-freedom of a field (component) defined on an
332: array of MOAB EntityHandles.
334: Not Collective
336: Input Parameters:
337: + dm - the discretization manager object
338: . npoints - the total number of Entities in the points array
339: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
340: - field - the field (component) index
342: Output Parameter:
343: . dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
345: Level: intermediate
347: .seealso: DMMoabGetFieldDof(), DMMoabGetFieldDofs()
348: @*/
349: PetscErrorCode DMMoabGetFieldDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt field, PetscInt* dof)
350: {
351: PetscInt i;
352: PetscErrorCode ierr;
353: DM_Moab *dmmoab;
358: dmmoab = (DM_Moab*)(dm)->data;
360: if (!dof) {
361: PetscMalloc1(npoints, &dof);
362: }
364: /* compute the DOF based on local blocking in the fields */
365: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
366: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
367: for (i = 0; i < npoints; ++i) {
368: dof[i] = (dmmoab->bs > 1 ? dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
369: dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + field * dmmoab->n);
370: }
371: return(0);
372: }
374: /*@C
375: DMMoabGetDofs - Gets the global degree-of-freedom for all fields (components) defined on an
376: array of MOAB EntityHandles.
378: Not Collective
380: Input Parameters:
381: + dm - the discretization manager object
382: . npoints - the total number of Entities in the points array
383: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
385: Output Parameter:
386: . dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
388: Level: intermediate
390: .seealso: DMMoabGetFieldDofs(), DMMoabGetDofsLocal(), DMMoabGetDofsBlocked()
391: @*/
392: PetscErrorCode DMMoabGetDofs(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
393: {
394: PetscInt i, field, offset;
395: PetscErrorCode ierr;
396: DM_Moab *dmmoab;
401: dmmoab = (DM_Moab*)(dm)->data;
403: if (!dof) {
404: PetscMalloc1(dmmoab->numFields * npoints, &dof);
405: }
407: /* compute the DOF based on local blocking in the fields */
408: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
409: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
410: for (field = 0; field < dmmoab->numFields; ++field) {
411: offset = field * dmmoab->n;
412: for (i = 0; i < npoints; ++i)
413: dof[i * dmmoab->numFields + field] = (dmmoab->bs > 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
414: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + offset);
415: }
416: return(0);
417: }
419: /*@C
420: DMMoabGetDofsLocal - Gets the local degree-of-freedom for all fields (components) defined on an
421: array of MOAB EntityHandles.
423: Not Collective
425: Input Parameters:
426: + dm - the discretization manager object
427: . npoints - the total number of Entities in the points array
428: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
430: Output Parameter:
431: . dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
433: Level: intermediate
435: .seealso: DMMoabGetFieldDofs(), DMMoabGetDofs(), DMMoabGetDofsBlocked()
436: @*/
437: PetscErrorCode DMMoabGetDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
438: {
439: PetscInt i, field, offset;
440: PetscErrorCode ierr;
441: DM_Moab *dmmoab;
446: dmmoab = (DM_Moab*)(dm)->data;
448: if (!dof) {
449: PetscMalloc1(dmmoab->numFields * npoints, &dof);
450: }
452: /* compute the DOF based on local blocking in the fields */
453: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
454: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
455: for (field = 0; field < dmmoab->numFields; ++field) {
456: offset = field * dmmoab->n;
457: for (i = 0; i < npoints; ++i)
458: dof[i * dmmoab->numFields + field] = (dmmoab->bs > 1 ? dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
459: dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + offset);
460: }
461: return(0);
462: }
464: /*@C
465: DMMoabGetDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
466: array of MOAB EntityHandles. It is useful when performing Blocked(Get/Set) methods in computation
467: of element residuals and assembly of the discrete systems when all fields are co-located.
469: Not Collective
471: Input Parameters:
472: + dm - the discretization manager object
473: . npoints - the total number of Entities in the points array
474: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
476: Output Parameter:
477: . dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat)
479: Level: intermediate
481: .seealso: DMMoabGetDofsLocal(), DMMoabGetDofs(), DMMoabGetDofsBlockedLocal()
482: @*/
483: PetscErrorCode DMMoabGetDofsBlocked(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
484: {
485: PetscInt i;
486: DM_Moab *dmmoab;
487: PetscErrorCode ierr;
492: dmmoab = (DM_Moab*)(dm)->data;
494: if (!dof) {
495: PetscMalloc1(npoints, &dof);
496: }
498: for (i = 0; i < npoints; ++i) {
499: dof[i] = dmmoab->gidmap[(PetscInt)points[i] - dmmoab->seqstart];
500: }
501: return(0);
502: }
504: /*@C
505: DMMoabGetDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
506: array of MOAB EntityHandles. It is useful when performing local Blocked(Get/Set) methods in computation
507: of element residuals and assembly of the discrete systems when all fields are co-located.
509: Not Collective
511: Input Parameters:
512: + dm - the discretization manager object
513: . npoints - the total number of Entities in the points array
514: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
516: Output Parameter:
517: . dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat)
519: Level: intermediate
521: .seealso: DMMoabGetDofsLocal(), DMMoabGetDofs(), DMMoabGetDofsBlockedLocal()
522: @*/
523: PetscErrorCode DMMoabGetDofsBlockedLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
524: {
525: PetscInt i;
526: DM_Moab *dmmoab;
527: PetscErrorCode ierr;
532: dmmoab = (DM_Moab*)(dm)->data;
534: if (!dof) {
535: PetscMalloc1(npoints, &dof);
536: }
538: for (i = 0; i < npoints; ++i)
539: dof[i] = dmmoab->lidmap[(PetscInt)points[i] - dmmoab->seqstart];
540: return(0);
541: }
543: /*@C
544: DMMoabGetVertexDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
545: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
546: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
548: Not Collective
550: Input Parameters:
551: . dm - the discretization manager object
553: Output Parameter:
554: . dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
556: Level: intermediate
558: .seealso: DMMoabGetVertexDofsBlockedLocal(), DMMoabGetDofsBlocked(), DMMoabGetDofsBlockedLocal()
559: @*/
560: PetscErrorCode DMMoabGetVertexDofsBlocked(DM dm, PetscInt** dof)
561: {
562: DM_Moab *dmmoab;
566: dmmoab = (DM_Moab*)(dm)->data;
568: *dof = dmmoab->gidmap;
569: return(0);
570: }
572: /*@C
573: DMMoabGetVertexDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
574: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
575: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
577: Not Collective
579: Input Parameters:
580: . dm - the discretization manager object
582: Output Parameter:
583: . dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
585: Level: intermediate
587: .seealso: DMMoabGetVertexDofsBlocked(), DMMoabGetDofsBlocked(), DMMoabGetDofsBlockedLocal()
588: @*/
589: PetscErrorCode DMMoabGetVertexDofsBlockedLocal(DM dm, PetscInt** dof)
590: {
591: DM_Moab *dmmoab;
596: dmmoab = (DM_Moab*)(dm)->data;
598: *dof = dmmoab->lidmap;
599: return(0);
600: }