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: .keywords: discretization manager, set, component solution
20: .seealso: DMMoabGetFieldName(), DMMoabSetGlobalFieldVector()
21: @*/
22: PetscErrorCodeDMMoabSetFieldVector(DM dm, PetscInt ifield, Vec fvec) 23: {
24: DM_Moab *dmmoab;
25: moab::Tag vtag, ntag;
26: const PetscScalar *varray;
27: PetscScalar *farray;
28: moab::ErrorCode merr;
29: PetscErrorCode ierr;
30: std::string tag_name;
34: dmmoab = (DM_Moab*)(dm)->data;
36: 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);
38: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
39: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag,
40: moab::MB_TAG_DENSE | moab::MB_TAG_CREAT); MBERRNM(merr);
42: DMMoabGetVecTag(fvec, &vtag);
44: merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
45: if (!tag_name.length() && merr != moab::MB_SUCCESS) {
46: VecGetArrayRead(fvec, &varray);
47: /* use the entity handle and the Dof index to set the right value */
48: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)varray); MBERRNM(merr);
49: VecRestoreArrayRead(fvec, &varray);
50: }
51: else {
52: PetscMalloc1(dmmoab->nloc, &farray);
53: /* we are using a MOAB Vec - directly copy the tag data to new one */
54: merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void*)farray); MBERRNM(merr);
55: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray); MBERRNM(merr);
56: /* make sure the parallel exchange for ghosts are done appropriately */
57: PetscFree(farray);
58: }
59: #ifdef MOAB_HAVE_MPI
60: merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vowned); MBERRNM(merr);
61: #endif
62: return(0);
63: }
66: /*@C
67: DMMoabSetGlobalFieldVector - Sets the vector reference that represents the global solution associated
68: with all fields (components) managed by DM.
69: A useful utility when updating the DM solution after a solve, to be serialized with the mesh for
70: checkpointing purposes.
72: Not Collective
74: Input Parameters:
75: + dm - the discretization manager object
76: . fvec - the global Vector solution corresponding to all the fields managed by DM 78: Level: intermediate
80: .keywords: discretization manager, set, component solution
82: .seealso: DMMoabGetFieldName(), DMMoabSetFieldVector()
83: @*/
84: PetscErrorCodeDMMoabSetGlobalFieldVector(DM dm, Vec fvec) 85: {
86: DM_Moab *dmmoab;
87: moab::Tag vtag, ntag;
88: const PetscScalar *rarray;
89: PetscScalar *varray, *farray;
90: moab::ErrorCode merr;
91: PetscErrorCode ierr;
92: PetscInt i, ifield;
93: std::string tag_name;
94: moab::Range::iterator iter;
98: dmmoab = (DM_Moab*)(dm)->data;
100: /* get the Tag corresponding to the global vector - possible that there is no tag associated.. */
101: DMMoabGetVecTag(fvec, &vtag);
102: merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
103: PetscMalloc1(dmmoab->nloc, &farray);
104: if (!tag_name.length() && merr != moab::MB_SUCCESS) {
105: /* not a MOAB vector - use VecGetSubVector to get the parts as needed */
106: VecGetArrayRead(fvec, &rarray);
107: for (ifield = 0; ifield < dmmoab->numFields; ++ifield) {
109: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
110: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag,
111: moab::MB_TAG_DENSE | moab::MB_TAG_CREAT); MBERRNM(merr);
113: for (i = 0; i < dmmoab->nloc; i++) {
114: farray[i] = (dmmoab->bs == 1 ? rarray[ifield * dmmoab->nloc + i] : rarray[i * dmmoab->numFields + ifield]);
115: }
117: /* use the entity handle and the Dof index to set the right value */
118: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray); MBERRNM(merr);
119: }
120: VecRestoreArrayRead(fvec, &rarray);
121: }
122: else {
123: PetscMalloc1(dmmoab->nloc * dmmoab->numFields, &varray);
125: /* we are using a MOAB Vec - directly copy the tag data to new one */
126: merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void*)varray); MBERRNM(merr);
127: for (ifield = 0; ifield < dmmoab->numFields; ++ifield) {
129: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
130: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag,
131: moab::MB_TAG_DENSE | moab::MB_TAG_CREAT); MBERRNM(merr);
133: /* we are using a MOAB Vec - directly copy the tag data to new one */
134: for (i = 0; i < dmmoab->nloc; i++) {
135: farray[i] = (dmmoab->bs == 1 ? varray[ifield * dmmoab->nloc + i] : varray[i * dmmoab->numFields + ifield]);
136: }
138: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray); MBERRNM(merr);
140: #ifdef MOAB_HAVE_MPI
141: /* make sure the parallel exchange for ghosts are done appropriately */
142: merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vlocal); MBERRNM(merr);
143: #endif
144: }
145: PetscFree(varray);
146: }
147: PetscFree(farray);
148: return(0);
149: }
152: /*@C
153: DMMoabSetFieldNames - Sets the number of fields and their names to be managed by the DM155: Not Collective
157: Input Parameters:
158: + dm - the discretization manager object
159: . numFields - the total number of fields
160: . fields - the array containing the names of each field (component); Can be NULL.
162: Level: intermediate
164: .keywords: discretization manager, set, component name
166: .seealso: DMMoabGetFieldName(), DMMoabSetFieldName()
167: @*/
168: PetscErrorCodeDMMoabSetFieldNames(DM dm, PetscInt numFields, const char* fields[])169: {
171: PetscInt i;
172: DM_Moab *dmmoab;
176: dmmoab = (DM_Moab*)(dm)->data;
178: /* first deallocate the existing field structure */
179: if (dmmoab->fieldNames) {
180: for (i = 0; i < dmmoab->numFields; i++) {
181: PetscFree(dmmoab->fieldNames[i]);
182: }
183: PetscFree(dmmoab->fieldNames);
184: }
186: /* now re-allocate and assign field names */
187: dmmoab->numFields = numFields;
188: PetscMalloc1(numFields, &dmmoab->fieldNames);
189: if (fields) {
190: for (i = 0; i < dmmoab->numFields; i++) {
191: PetscStrallocpy(fields[i], (char**) &dmmoab->fieldNames[i]);
192: }
193: }
194: DMSetNumFields(dm, numFields);
195: return(0);
196: }
199: /*@C
200: DMMoabGetFieldName - Gets the names of individual field components in multicomponent
201: vectors associated with a DMDA.
203: Not Collective
205: Input Parameter:
206: + dm - the discretization manager object
207: . field - field number for the DMMoab (0, 1, ... dof-1), where dof indicates the
208: number of degrees of freedom per node within the DMMoab
210: Output Parameter:
211: . fieldName - the name of the field (component)
213: Level: intermediate
215: .keywords: discretization manager, get, component name
217: .seealso: DMMoabSetFieldName(), DMMoabSetFields()
218: @*/
219: PetscErrorCodeDMMoabGetFieldName(DM dm, PetscInt field, const char **fieldName)220: {
221: DM_Moab *dmmoab;
225: dmmoab = (DM_Moab*)(dm)->data;
226: 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);
228: *fieldName = dmmoab->fieldNames[field];
229: return(0);
230: }
233: /*@C
234: DMMoabSetFieldName - Sets the name of a field (component) managed by the DM236: Not Collective
238: Input Parameters:
239: + dm - the discretization manager object
240: . field - the field number
241: . fieldName - the field (component) name
243: Level: intermediate
244: Notes: Can only be called after DMMoabSetFields supplied with correct numFields
246: .keywords: discretization manager, set, component name
248: .seealso: DMMoabGetFieldName(), DMMoabSetFields()
249: @*/
250: PetscErrorCodeDMMoabSetFieldName(DM dm, PetscInt field, const char *fieldName)251: {
253: DM_Moab *dmmoab;
259: dmmoab = (DM_Moab*)(dm)->data;
260: 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);
262: if (dmmoab->fieldNames[field]) {
263: PetscFree(dmmoab->fieldNames[field]);
264: }
265: PetscStrallocpy(fieldName, (char**) &dmmoab->fieldNames[field]);
266: return(0);
267: }
270: /*@C
271: DMMoabGetFieldDof - Gets the global degree-of-freedom of a field (component) defined on a
272: particular MOAB EntityHandle.
274: Not Collective
276: Input Parameters:
277: + dm - the discretization manager object
278: . point - the MOAB EntityHandle container which holds the field degree-of-freedom values
279: . field - the field (component) index
281: Output Parameter:
282: + dof - the global degree-of-freedom index corresponding to the field in the discrete representation (Vec, Mat)
284: Level: beginner
286: .keywords: discretization manager, get, global degree of freedom
288: .seealso: DMMoabGetFieldDofs(), DMMoabGetFieldDofsLocal()
289: @*/
290: PetscErrorCodeDMMoabGetFieldDof(DM dm, moab::EntityHandle point, PetscInt field, PetscInt* dof)
291: {
292: DM_Moab *dmmoab;
296: dmmoab = (DM_Moab*)(dm)->data;
298: *dof = (dmmoab->bs == 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(point) - dmmoab->seqstart] + field * dmmoab->n :
299: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(point) - dmmoab->seqstart] * dmmoab->numFields + field);
300: return(0);
301: }
304: /*@C
305: DMMoabGetFieldDofs - Gets the global degree-of-freedom of a field (component) defined on an
306: array of MOAB EntityHandles.
308: Not Collective
310: Input Parameters:
311: + dm - the discretization manager object
312: . npoints - the total number of Entities in the points array
313: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
314: . field - the field (component) index
316: Output Parameter:
317: + dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
319: Level: intermediate
321: .keywords: discretization manager, get, global degrees of freedom
323: .seealso: DMMoabGetFieldDof(), DMMoabGetFieldDofsLocal()
324: @*/
325: PetscErrorCodeDMMoabGetFieldDofs(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt field, PetscInt* dof)
326: {
327: PetscInt i;
328: PetscErrorCode ierr;
329: DM_Moab *dmmoab;
334: dmmoab = (DM_Moab*)(dm)->data;
336: if (!dof) {
337: PetscMalloc1(npoints, &dof);
338: }
340: /* compute the DOF based on local blocking in the fields */
341: /* We also assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
342: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
343: for (i = 0; i < npoints; ++i)
344: dof[i] = (dmmoab->bs == 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + field * dmmoab->n :
345: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field);
346: return(0);
347: }
350: /*@C
351: DMMoabGetFieldDofsLocal - Gets the local degrees-of-freedom of a field (component) defined on an
352: array of MOAB EntityHandles.
354: Not Collective
356: Input Parameters:
357: + dm - the discretization manager object
358: . npoints - the total number of Entities in the points array
359: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
360: . field - the field (component) index
362: Output Parameter:
363: + dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
365: Level: intermediate
367: .keywords: discretization manager, get, local degrees of freedom
369: .seealso: DMMoabGetFieldDof(), DMMoabGetFieldDofs()
370: @*/
371: PetscErrorCodeDMMoabGetFieldDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt field, PetscInt* dof)
372: {
373: PetscInt i;
374: PetscErrorCode ierr;
375: DM_Moab *dmmoab;
380: dmmoab = (DM_Moab*)(dm)->data;
382: if (!dof) {
383: PetscMalloc1(npoints, &dof);
384: }
386: /* compute the DOF based on local blocking in the fields */
387: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
388: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
389: for (i = 0; i < npoints; ++i) {
390: dof[i] = (dmmoab->bs > 1 ? dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
391: dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + field * dmmoab->n);
392: }
393: return(0);
394: }
397: /*@C
398: DMMoabGetDofs - Gets the global degree-of-freedom for all fields (components) defined on an
399: array of MOAB EntityHandles.
401: Not Collective
403: Input Parameters:
404: + dm - the discretization manager object
405: . npoints - the total number of Entities in the points array
406: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
408: Output Parameter:
409: + dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
411: Level: intermediate
413: .keywords: discretization manager, get, global degrees of freedom
415: .seealso: DMMoabGetFieldDofs(), DMMoabGetDofsLocal(), DMMoabGetDofsBlocked()
416: @*/
417: PetscErrorCodeDMMoabGetDofs(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
418: {
419: PetscInt i, field, offset;
420: PetscErrorCode ierr;
421: DM_Moab *dmmoab;
426: dmmoab = (DM_Moab*)(dm)->data;
428: if (!dof) {
429: PetscMalloc1(dmmoab->numFields * npoints, &dof);
430: }
432: /* compute the DOF based on local blocking in the fields */
433: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
434: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
435: for (field = 0; field < dmmoab->numFields; ++field) {
436: offset = field * dmmoab->n;
437: for (i = 0; i < npoints; ++i)
438: dof[i * dmmoab->numFields + field] = (dmmoab->bs > 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
439: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + offset);
440: }
441: return(0);
442: }
444: /*@C
445: DMMoabGetDofsLocal - Gets the local degree-of-freedom for all fields (components) defined on an
446: array of MOAB EntityHandles.
448: Not Collective
450: Input Parameters:
451: + dm - the discretization manager object
452: . npoints - the total number of Entities in the points array
453: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
455: Output Parameter:
456: + dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
458: Level: intermediate
460: .keywords: discretization manager, get, global degrees of freedom
462: .seealso: DMMoabGetFieldDofs(), DMMoabGetDofs(), DMMoabGetDofsBlocked()
463: @*/
464: PetscErrorCodeDMMoabGetDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
465: {
466: PetscInt i, field, offset;
467: PetscErrorCode ierr;
468: DM_Moab *dmmoab;
473: dmmoab = (DM_Moab*)(dm)->data;
475: if (!dof) {
476: PetscMalloc1(dmmoab->numFields * npoints, &dof);
477: }
479: /* compute the DOF based on local blocking in the fields */
480: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
481: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
482: for (field = 0; field < dmmoab->numFields; ++field) {
483: offset = field * dmmoab->n;
484: for (i = 0; i < npoints; ++i)
485: dof[i * dmmoab->numFields + field] = (dmmoab->bs > 1 ? dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
486: dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + offset);
487: }
488: return(0);
489: }
492: /*@C
493: DMMoabGetDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
494: array of MOAB EntityHandles. It is useful when performing Blocked(Get/Set) methods in computation
495: of element residuals and assembly of the discrete systems when all fields are co-located.
497: Not Collective
499: Input Parameters:
500: + dm - the discretization manager object
501: . npoints - the total number of Entities in the points array
502: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
504: Output Parameter:
505: + dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat)
507: Level: intermediate
509: .keywords: discretization manager, get, global degrees of freedom
511: .seealso: DMMoabGetDofsLocal(), DMMoabGetDofs(), DMMoabGetDofsBlockedLocal()
512: @*/
513: PetscErrorCodeDMMoabGetDofsBlocked(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
514: {
515: PetscInt i;
516: DM_Moab *dmmoab;
517: PetscErrorCode ierr;
522: dmmoab = (DM_Moab*)(dm)->data;
524: if (!dof) {
525: PetscMalloc1(npoints, &dof);
526: }
528: for (i = 0; i < npoints; ++i) {
529: dof[i] = dmmoab->gidmap[(PetscInt)points[i] - dmmoab->seqstart];
530: }
531: return(0);
532: }
535: /*@C
536: DMMoabGetDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
537: array of MOAB EntityHandles. It is useful when performing local Blocked(Get/Set) methods in computation
538: of element residuals and assembly of the discrete systems when all fields are co-located.
540: Not Collective
542: Input Parameters:
543: + dm - the discretization manager object
544: . npoints - the total number of Entities in the points array
545: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
547: Output Parameter:
548: + dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat)
550: Level: intermediate
552: .keywords: discretization manager, get, global degrees of freedom
554: .seealso: DMMoabGetDofsLocal(), DMMoabGetDofs(), DMMoabGetDofsBlockedLocal()
555: @*/
556: PetscErrorCodeDMMoabGetDofsBlockedLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
557: {
558: PetscInt i;
559: DM_Moab *dmmoab;
560: PetscErrorCode ierr;
565: dmmoab = (DM_Moab*)(dm)->data;
567: if (!dof) {
568: PetscMalloc1(npoints, &dof);
569: }
571: for (i = 0; i < npoints; ++i)
572: dof[i] = dmmoab->lidmap[(PetscInt)points[i] - dmmoab->seqstart];
573: return(0);
574: }
577: /*@C
578: DMMoabGetVertexDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
579: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
580: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
582: Not Collective
584: Input Parameters:
585: + dm - the discretization manager object
587: Output Parameter:
588: + dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
590: Level: intermediate
592: .keywords: discretization manager, get, blocked degrees of freedom
594: .seealso: DMMoabGetVertexDofsBlockedLocal(), DMMoabGetDofsBlocked(), DMMoabGetDofsBlockedLocal()
595: @*/
596: PetscErrorCodeDMMoabGetVertexDofsBlocked(DM dm, PetscInt** dof)597: {
598: DM_Moab *dmmoab;
602: dmmoab = (DM_Moab*)(dm)->data;
604: *dof = dmmoab->gidmap;
605: return(0);
606: }
609: /*@C
610: DMMoabGetVertexDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
611: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
612: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
614: Not Collective
616: Input Parameters:
617: + dm - the discretization manager object
619: Output Parameter:
620: + dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
622: Level: intermediate
624: .keywords: discretization manager, get, blocked degrees of freedom
626: .seealso: DMMoabGetVertexDofsBlocked(), DMMoabGetDofsBlocked(), DMMoabGetDofsBlockedLocal()
627: @*/
628: PetscErrorCodeDMMoabGetVertexDofsBlockedLocal(DM dm, PetscInt** dof)629: {
630: DM_Moab *dmmoab;
635: dmmoab = (DM_Moab*)(dm)->data;
637: *dof = dmmoab->lidmap;
638: return(0);
639: }