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:
245: Can only be called after DMMoabSetFields supplied with correct numFields
247: .keywords: discretization manager, set, component name
249: .seealso: DMMoabGetFieldName(), DMMoabSetFields()
250: @*/
251: PetscErrorCodeDMMoabSetFieldName(DM dm, PetscInt field, const char *fieldName)252: {
254: DM_Moab *dmmoab;
260: dmmoab = (DM_Moab*)(dm)->data;
261: 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);
263: if (dmmoab->fieldNames[field]) {
264: PetscFree(dmmoab->fieldNames[field]);
265: }
266: PetscStrallocpy(fieldName, (char**) &dmmoab->fieldNames[field]);
267: return(0);
268: }
271: /*@C
272: DMMoabGetFieldDof - Gets the global degree-of-freedom of a field (component) defined on a
273: particular MOAB EntityHandle.
275: Not Collective
277: Input Parameters:
278: + dm - the discretization manager object
279: . point - the MOAB EntityHandle container 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 corresponding to the field in the discrete representation (Vec, Mat)
285: Level: beginner
287: .keywords: discretization manager, get, global degree of freedom
289: .seealso: DMMoabGetFieldDofs(), DMMoabGetFieldDofsLocal()
290: @*/
291: PetscErrorCodeDMMoabGetFieldDof(DM dm, moab::EntityHandle point, PetscInt field, PetscInt* dof)
292: {
293: DM_Moab *dmmoab;
297: dmmoab = (DM_Moab*)(dm)->data;
299: *dof = (dmmoab->bs == 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(point) - dmmoab->seqstart] + field * dmmoab->n :
300: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(point) - dmmoab->seqstart] * dmmoab->numFields + field);
301: return(0);
302: }
305: /*@C
306: DMMoabGetFieldDofs - Gets the global degree-of-freedom of a field (component) defined on an
307: array of MOAB EntityHandles.
309: Not Collective
311: Input Parameters:
312: + dm - the discretization manager object
313: . npoints - the total number of Entities in the points array
314: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
315: . field - the field (component) index
317: Output Parameter:
318: + dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
320: Level: intermediate
322: .keywords: discretization manager, get, global degrees of freedom
324: .seealso: DMMoabGetFieldDof(), DMMoabGetFieldDofsLocal()
325: @*/
326: PetscErrorCodeDMMoabGetFieldDofs(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt field, PetscInt* dof)
327: {
328: PetscInt i;
329: PetscErrorCode ierr;
330: DM_Moab *dmmoab;
335: dmmoab = (DM_Moab*)(dm)->data;
337: if (!dof) {
338: PetscMalloc1(npoints, &dof);
339: }
341: /* compute the DOF based on local blocking in the fields */
342: /* We also assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
343: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
344: for (i = 0; i < npoints; ++i)
345: dof[i] = (dmmoab->bs == 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + field * dmmoab->n :
346: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field);
347: return(0);
348: }
351: /*@C
352: DMMoabGetFieldDofsLocal - Gets the local degrees-of-freedom of a field (component) defined on an
353: array of MOAB EntityHandles.
355: Not Collective
357: Input Parameters:
358: + dm - the discretization manager object
359: . npoints - the total number of Entities in the points array
360: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
361: . field - the field (component) index
363: Output Parameter:
364: + dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
366: Level: intermediate
368: .keywords: discretization manager, get, local degrees of freedom
370: .seealso: DMMoabGetFieldDof(), DMMoabGetFieldDofs()
371: @*/
372: PetscErrorCodeDMMoabGetFieldDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt field, PetscInt* dof)
373: {
374: PetscInt i;
375: PetscErrorCode ierr;
376: DM_Moab *dmmoab;
381: dmmoab = (DM_Moab*)(dm)->data;
383: if (!dof) {
384: PetscMalloc1(npoints, &dof);
385: }
387: /* compute the DOF based on local blocking in the fields */
388: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
389: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
390: for (i = 0; i < npoints; ++i) {
391: dof[i] = (dmmoab->bs > 1 ? dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
392: dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + field * dmmoab->n);
393: }
394: return(0);
395: }
398: /*@C
399: DMMoabGetDofs - Gets the global degree-of-freedom for all fields (components) defined on an
400: array of MOAB EntityHandles.
402: Not Collective
404: Input Parameters:
405: + dm - the discretization manager object
406: . npoints - the total number of Entities in the points array
407: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
409: Output Parameter:
410: + dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
412: Level: intermediate
414: .keywords: discretization manager, get, global degrees of freedom
416: .seealso: DMMoabGetFieldDofs(), DMMoabGetDofsLocal(), DMMoabGetDofsBlocked()
417: @*/
418: PetscErrorCodeDMMoabGetDofs(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
419: {
420: PetscInt i, field, offset;
421: PetscErrorCode ierr;
422: DM_Moab *dmmoab;
427: dmmoab = (DM_Moab*)(dm)->data;
429: if (!dof) {
430: PetscMalloc1(dmmoab->numFields * npoints, &dof);
431: }
433: /* compute the DOF based on local blocking in the fields */
434: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
435: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
436: for (field = 0; field < dmmoab->numFields; ++field) {
437: offset = field * dmmoab->n;
438: for (i = 0; i < npoints; ++i)
439: dof[i * dmmoab->numFields + field] = (dmmoab->bs > 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
440: dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + offset);
441: }
442: return(0);
443: }
445: /*@C
446: DMMoabGetDofsLocal - Gets the local degree-of-freedom for all fields (components) defined on an
447: array of MOAB EntityHandles.
449: Not Collective
451: Input Parameters:
452: + dm - the discretization manager object
453: . npoints - the total number of Entities in the points array
454: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
456: Output Parameter:
457: + dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
459: Level: intermediate
461: .keywords: discretization manager, get, global degrees of freedom
463: .seealso: DMMoabGetFieldDofs(), DMMoabGetDofs(), DMMoabGetDofsBlocked()
464: @*/
465: PetscErrorCodeDMMoabGetDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
466: {
467: PetscInt i, field, offset;
468: PetscErrorCode ierr;
469: DM_Moab *dmmoab;
474: dmmoab = (DM_Moab*)(dm)->data;
476: if (!dof) {
477: PetscMalloc1(dmmoab->numFields * npoints, &dof);
478: }
480: /* compute the DOF based on local blocking in the fields */
481: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
482: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
483: for (field = 0; field < dmmoab->numFields; ++field) {
484: offset = field * dmmoab->n;
485: for (i = 0; i < npoints; ++i)
486: dof[i * dmmoab->numFields + field] = (dmmoab->bs > 1 ? dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field :
487: dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + offset);
488: }
489: return(0);
490: }
493: /*@C
494: DMMoabGetDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
495: array of MOAB EntityHandles. It is useful when performing Blocked(Get/Set) methods in computation
496: of element residuals and assembly of the discrete systems when all fields are co-located.
498: Not Collective
500: Input Parameters:
501: + dm - the discretization manager object
502: . npoints - the total number of Entities in the points array
503: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
505: Output Parameter:
506: + dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat)
508: Level: intermediate
510: .keywords: discretization manager, get, global degrees of freedom
512: .seealso: DMMoabGetDofsLocal(), DMMoabGetDofs(), DMMoabGetDofsBlockedLocal()
513: @*/
514: PetscErrorCodeDMMoabGetDofsBlocked(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
515: {
516: PetscInt i;
517: DM_Moab *dmmoab;
518: PetscErrorCode ierr;
523: dmmoab = (DM_Moab*)(dm)->data;
525: if (!dof) {
526: PetscMalloc1(npoints, &dof);
527: }
529: for (i = 0; i < npoints; ++i) {
530: dof[i] = dmmoab->gidmap[(PetscInt)points[i] - dmmoab->seqstart];
531: }
532: return(0);
533: }
536: /*@C
537: DMMoabGetDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
538: array of MOAB EntityHandles. It is useful when performing local Blocked(Get/Set) methods in computation
539: of element residuals and assembly of the discrete systems when all fields are co-located.
541: Not Collective
543: Input Parameters:
544: + dm - the discretization manager object
545: . npoints - the total number of Entities in the points array
546: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
548: Output Parameter:
549: + dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat)
551: Level: intermediate
553: .keywords: discretization manager, get, global degrees of freedom
555: .seealso: DMMoabGetDofsLocal(), DMMoabGetDofs(), DMMoabGetDofsBlockedLocal()
556: @*/
557: PetscErrorCodeDMMoabGetDofsBlockedLocal(DM dm, PetscInt npoints, const moab::EntityHandle* points, PetscInt* dof)
558: {
559: PetscInt i;
560: DM_Moab *dmmoab;
561: PetscErrorCode ierr;
566: dmmoab = (DM_Moab*)(dm)->data;
568: if (!dof) {
569: PetscMalloc1(npoints, &dof);
570: }
572: for (i = 0; i < npoints; ++i)
573: dof[i] = dmmoab->lidmap[(PetscInt)points[i] - dmmoab->seqstart];
574: return(0);
575: }
578: /*@C
579: DMMoabGetVertexDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
580: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
581: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
583: Not Collective
585: Input Parameters:
586: + dm - the discretization manager object
588: Output Parameter:
589: + dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
591: Level: intermediate
593: .keywords: discretization manager, get, blocked degrees of freedom
595: .seealso: DMMoabGetVertexDofsBlockedLocal(), DMMoabGetDofsBlocked(), DMMoabGetDofsBlockedLocal()
596: @*/
597: PetscErrorCodeDMMoabGetVertexDofsBlocked(DM dm, PetscInt** dof)598: {
599: DM_Moab *dmmoab;
603: dmmoab = (DM_Moab*)(dm)->data;
605: *dof = dmmoab->gidmap;
606: return(0);
607: }
610: /*@C
611: DMMoabGetVertexDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
612: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
613: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
615: Not Collective
617: Input Parameters:
618: + dm - the discretization manager object
620: Output Parameter:
621: + dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
623: Level: intermediate
625: .keywords: discretization manager, get, blocked degrees of freedom
627: .seealso: DMMoabGetVertexDofsBlocked(), DMMoabGetDofsBlocked(), DMMoabGetDofsBlockedLocal()
628: @*/
629: PetscErrorCodeDMMoabGetVertexDofsBlockedLocal(DM dm, PetscInt** dof)630: {
631: DM_Moab *dmmoab;
636: dmmoab = (DM_Moab*)(dm)->data;
638: *dof = dmmoab->lidmap;
639: return(0);
640: }