1: #include <petsc/private/dmmbimpl.h> /*I "petscdmmoab.h" I*/
3: #include <petscdmmoab.h>
7: /*@C
8: DMMoabSetFieldVector - Sets the vector reference that represents the solution associated
9: with a particular field component.
11: Not Collective
13: Input Parameters:
14: + dm - the discretization manager object
15: . ifield - the index of the field as set before via DMMoabSetFieldName.
16: . fvec - the Vector solution corresponding to the field (component)
18: Level: intermediate
20: .keywords: discretization manager, set, component solution
22: .seealso: DMMoabGetFieldName(), DMMoabSetGlobalFieldVector()
23: @*/
24: PetscErrorCodeDMMoabSetFieldVector(DM dm, PetscInt ifield, Vec fvec) 25: {
26: DM_Moab *dmmoab;
27: moab::Tag vtag,ntag;
28: const PetscScalar *varray;
29: PetscScalar *farray;
30: moab::ErrorCode merr;
31: PetscErrorCode ierr;
32: std::string tag_name;
36: dmmoab = (DM_Moab*)(dm)->data;
38: 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);
40: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
41: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield],1,moab::MB_TYPE_DOUBLE,ntag,
42: moab::MB_TAG_DENSE|moab::MB_TAG_CREAT);MBERRNM(merr);
44: DMMoabGetVecTag(fvec,&vtag);
46: merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
47: if (!tag_name.length() && merr !=moab::MB_SUCCESS) {
48: VecGetArrayRead(fvec,&varray);
49: /* use the entity handle and the Dof index to set the right value */
50: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)varray);MBERRNM(merr);
51: VecRestoreArrayRead(fvec,&varray);
52: }
53: else {
54: PetscMalloc1(dmmoab->nloc,&farray);
55: /* we are using a MOAB Vec - directly copy the tag data to new one */
56: merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void*)farray);MBERRNM(merr);
57: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray);MBERRNM(merr);
58: /* make sure the parallel exchange for ghosts are done appropriately */
59: PetscFree(farray);
60: }
61: merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vowned);MBERRNM(merr);
62: return(0);
63: }
68: /*@C
69: DMMoabSetGlobalFieldVector - Sets the vector reference that represents the global solution associated
70: with all fields (components) managed by DM.
71: A useful utility when updating the DM solution after a solve, to be serialized with the mesh for
72: checkpointing purposes.
74: Not Collective
76: Input Parameters:
77: + dm - the discretization manager object
78: . fvec - the global Vector solution corresponding to all the fields managed by DM 80: Level: intermediate
82: .keywords: discretization manager, set, component solution
84: .seealso: DMMoabGetFieldName(), DMMoabSetFieldVector()
85: @*/
86: PetscErrorCodeDMMoabSetGlobalFieldVector(DM dm, Vec fvec) 87: {
88: DM_Moab *dmmoab;
89: moab::Tag vtag,ntag;
90: const PetscScalar *rarray;
91: PetscScalar *varray,*farray;
92: moab::ErrorCode merr;
93: PetscErrorCode ierr;
94: PetscInt i,ifield;
95: std::string tag_name;
96: moab::Range::iterator iter;
100: dmmoab = (DM_Moab*)(dm)->data;
102: /* get the Tag corresponding to the global vector - possible that there is no tag associated.. */
103: DMMoabGetVecTag(fvec,&vtag);
104: merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
105: PetscMalloc1(dmmoab->nloc,&farray);
106: if (!tag_name.length() && merr !=moab::MB_SUCCESS) {
107: /* not a MOAB vector - use VecGetSubVector to get the parts as needed */
109: VecGetArrayRead(fvec,&rarray);
110: for (ifield=0; ifield<dmmoab->numFields; ++ifield) {
112: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
113: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield],1,moab::MB_TYPE_DOUBLE,ntag,
114: moab::MB_TAG_DENSE|moab::MB_TAG_CREAT);MBERRNM(merr);
116: for(i=0;i<dmmoab->nloc;i++) {
117: if (dmmoab->bs == 1)
118: farray[i]=rarray[ifield*dmmoab->nloc+i];
119: else
120: farray[i]=rarray[i*dmmoab->numFields+ifield];
121: }
123: /* use the entity handle and the Dof index to set the right value */
124: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray);MBERRNM(merr);
125: }
126: VecRestoreArrayRead(fvec,&rarray);
127: }
128: else {
129: PetscMalloc1(dmmoab->nloc*dmmoab->numFields,&varray);
131: /* we are using a MOAB Vec - directly copy the tag data to new one */
132: merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void*)varray);MBERRNM(merr);
133: for (ifield=0; ifield<dmmoab->numFields; ++ifield) {
135: /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
136: merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield],1,moab::MB_TYPE_DOUBLE,ntag,
137: moab::MB_TAG_DENSE|moab::MB_TAG_CREAT);MBERRNM(merr);
139: /* we are using a MOAB Vec - directly copy the tag data to new one */
140: for(i=0; i < dmmoab->nloc; i++) {
141: if (dmmoab->bs == 1)
142: farray[i]=varray[ifield*dmmoab->nloc+i];
143: else
144: farray[i]=varray[i*dmmoab->numFields+ifield];
145: }
147: merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void*)farray);MBERRNM(merr);
148: /* make sure the parallel exchange for ghosts are done appropriately */
149: merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vlocal);MBERRNM(merr);
150: }
151: PetscFree(varray);
152: }
153: PetscFree(farray);
154: return(0);
155: }
160: /*@C
161: DMMoabSetFieldNames - Sets the number of fields and their names to be managed by the DM163: Not Collective
165: Input Parameters:
166: + dm - the discretization manager object
167: . numFields - the total number of fields
168: . fields - the array containing the names of each field (component); Can be NULL.
170: Level: intermediate
172: .keywords: discretization manager, set, component name
173: 174: .seealso: DMMoabGetFieldName(), DMMoabSetFieldName()
175: @*/
176: PetscErrorCodeDMMoabSetFieldNames(DM dm,PetscInt numFields,const char* fields[])177: {
179: PetscInt i;
180: DM_Moab *dmmoab;
184: dmmoab = (DM_Moab*)(dm)->data;
186: /* first deallocate the existing field structure */
187: if (dmmoab->fieldNames) {
188: for(i=0; i<dmmoab->numFields; i++) {
189: PetscFree(dmmoab->fieldNames[i]);
190: }
191: PetscFree(dmmoab->fieldNames);
192: }
194: /* now re-allocate and assign field names */
195: dmmoab->numFields = numFields;
196: PetscMalloc(numFields*sizeof(char*),&dmmoab->fieldNames);
197: if (fields) {
198: for(i=0; i<dmmoab->numFields; i++) {
199: PetscStrallocpy(fields[i], (char**) &dmmoab->fieldNames[i]);
200: }
201: }
202: return(0);
203: }
208: /*@C
209: DMMoabGetFieldName - Gets the names of individual field components in multicomponent
210: vectors associated with a DMDA.
212: Not Collective
214: Input Parameter:
215: + dm - the discretization manager object
216: . field - field number for the DMMoab (0, 1, ... dof-1), where dof indicates the
217: number of degrees of freedom per node within the DMMoab
219: Output Parameter:
220: . fieldName - the name of the field (component)
222: Level: intermediate
224: .keywords: discretization manager, get, component name
226: .seealso: DMMoabSetFieldName(), DMMoabSetFields()
227: @*/
228: PetscErrorCodeDMMoabGetFieldName(DM dm, PetscInt field, const char **fieldName)229: {
230: DM_Moab *dmmoab;
234: dmmoab = (DM_Moab*)(dm)->data;
235: 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);
237: *fieldName = dmmoab->fieldNames[field];
238: return(0);
239: }
244: /*@C
245: DMMoabSetFieldName - Sets the name of a field (component) managed by the DM247: Not Collective
249: Input Parameters:
250: + dm - the discretization manager object
251: . field - the field number
252: . fieldName - the field (component) name
254: Level: intermediate
255: Notes: Can only be called after DMMoabSetFields supplied with correct numFields
257: .keywords: discretization manager, set, component name
258: 259: .seealso: DMMoabGetFieldName(), DMMoabSetFields()
260: @*/
261: PetscErrorCodeDMMoabSetFieldName(DM dm, PetscInt field, const char *fieldName)262: {
264: DM_Moab *dmmoab;
270: dmmoab = (DM_Moab*)(dm)->data;
271: 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);
273: if (dmmoab->fieldNames[field]) {
274: PetscFree(dmmoab->fieldNames[field]);
275: }
276: PetscStrallocpy(fieldName, (char**) &dmmoab->fieldNames[field]);
277: return(0);
278: }
283: /*@C
284: DMMoabGetFieldDof - Gets the global degree-of-freedom of a field (component) defined on a
285: particular MOAB EntityHandle.
287: Not Collective
289: Input Parameters:
290: + dm - the discretization manager object
291: . point - the MOAB EntityHandle container which holds the field degree-of-freedom values
292: . field - the field (component) index
294: Output Parameter:
295: + dof - the global degree-of-freedom index corresponding to the field in the discrete representation (Vec, Mat)
297: Level: beginner
299: .keywords: discretization manager, get, global degree of freedom
301: .seealso: DMMoabGetFieldDofs(), DMMoabGetFieldDofsLocal()
302: @*/
303: PetscErrorCodeDMMoabGetFieldDof(DM dm,moab::EntityHandle point,PetscInt field,PetscInt* dof)
304: {
305: DM_Moab *dmmoab;
309: dmmoab = (DM_Moab*)(dm)->data;
311: *dof=dmmoab->gidmap[(PetscInt)point];
312: return(0);
313: }
318: /*@C
319: DMMoabGetFieldDofs - Gets the global degree-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 global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
333: Level: intermediate
335: .keywords: discretization manager, get, global degrees of freedom
336: 337: .seealso: DMMoabGetFieldDof(), DMMoabGetFieldDofsLocal()
338: @*/
339: PetscErrorCodeDMMoabGetFieldDofs(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt field,PetscInt* dof)
340: {
341: PetscInt i;
342: PetscErrorCode ierr;
343: DM_Moab *dmmoab;
348: dmmoab = (DM_Moab*)(dm)->data;
350: if (!dof) {
351: PetscMalloc1(npoints,&dof);
352: }
354: if (dmmoab->bs > 1) { /* compute the DOF based on local blocking in the fields */
355: for (i=0; i<npoints; ++i)
356: dof[i] = dmmoab->gidmap[(PetscInt)points[i]]*dmmoab->numFields+field;
357: }
358: else {
359: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
360: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
361: for (i=0; i<npoints; ++i)
362: dof[i] = dmmoab->gidmap[(PetscInt)points[i]]+field*dmmoab->n;
363: }
364: return(0);
365: }
370: /*@C
371: DMMoabGetFieldDofsLocal - Gets the local degrees-of-freedom of a field (component) defined on an
372: array of MOAB EntityHandles.
374: Not Collective
376: Input Parameters:
377: + dm - the discretization manager object
378: . npoints - the total number of Entities in the points array
379: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
380: . field - the field (component) index
382: Output Parameter:
383: + dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
385: Level: intermediate
387: .keywords: discretization manager, get, local degrees of freedom
388: 389: .seealso: DMMoabGetFieldDof(), DMMoabGetFieldDofs()
390: @*/
391: PetscErrorCodeDMMoabGetFieldDofsLocal(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt field,PetscInt* dof)
392: {
393: PetscInt i,offset;
394: PetscErrorCode ierr;
395: DM_Moab *dmmoab;
400: dmmoab = (DM_Moab*)(dm)->data;
402: if (!dof) {
403: PetscMalloc1(npoints,&dof);
404: }
406: if (dmmoab->bs > 1) { /* compute the DOF based on local blocking in the fields */
407: for (i=0; i<npoints; ++i)
408: dof[i] = dmmoab->lidmap[(PetscInt)points[i]]*dmmoab->numFields+field;
409: }
410: else {
411: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
412: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
413: offset = field*dmmoab->n;
414: for (i=0; i<npoints; ++i)
415: dof[i] = dmmoab->lidmap[(PetscInt)points[i]]+offset;
416: }
417: return(0);
418: }
423: /*@C
424: DMMoabGetDofs - Gets the global degree-of-freedom for all fields (components) defined on an
425: array of MOAB EntityHandles.
427: Not Collective
429: Input Parameters:
430: + dm - the discretization manager object
431: . npoints - the total number of Entities in the points array
432: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
434: Output Parameter:
435: + dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
437: Level: intermediate
439: .keywords: discretization manager, get, global degrees of freedom
440: 441: .seealso: DMMoabGetFieldDofs(), DMMoabGetDofsLocal(), DMMoabGetDofsBlocked()
442: @*/
443: PetscErrorCodeDMMoabGetDofs(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof)
444: {
445: PetscInt i,field,offset;
446: PetscErrorCode ierr;
447: DM_Moab *dmmoab;
452: dmmoab = (DM_Moab*)(dm)->data;
454: if (!dof) {
455: PetscMalloc1(dmmoab->numFields*npoints,&dof);
456: }
457: 458: if (dmmoab->bs > 1) { /* compute the DOF based on local blocking in the fields */
459: for (field=0; field<dmmoab->numFields; ++field) {
460: for (i=0; i<npoints; ++i)
461: dof[i*dmmoab->numFields+field] = dmmoab->gidmap[(PetscInt)points[i]]*dmmoab->numFields+field;
462: }
463: }
464: else {
465: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
466: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
467: for (field=0; field<dmmoab->numFields; ++field) {
468: offset = field*dmmoab->n;
469: for (i=0; i<npoints; ++i)
470: dof[i*dmmoab->numFields+field] = dmmoab->gidmap[(PetscInt)points[i]]+offset;
471: }
472: }
473: return(0);
474: }
479: /*@C
480: DMMoabGetDofsLocal - Gets the local degree-of-freedom for all fields (components) defined on an
481: array of MOAB EntityHandles.
483: Not Collective
485: Input Parameters:
486: + dm - the discretization manager object
487: . npoints - the total number of Entities in the points array
488: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
490: Output Parameter:
491: + dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)
493: Level: intermediate
495: .keywords: discretization manager, get, global degrees of freedom
496: 497: .seealso: DMMoabGetFieldDofs(), DMMoabGetDofs(), DMMoabGetDofsBlocked()
498: @*/
499: PetscErrorCodeDMMoabGetDofsLocal(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof)
500: {
501: PetscInt i,field,offset;
502: PetscErrorCode ierr;
503: DM_Moab *dmmoab;
508: dmmoab = (DM_Moab*)(dm)->data;
510: if (!dof) {
511: PetscMalloc1(dmmoab->numFields*npoints,&dof);
512: }
514: if (dmmoab->bs > 1) { /* compute the DOF based on local blocking in the fields */
515: for (field=0; field<dmmoab->numFields; ++field) {
516: for (i=0; i<npoints; ++i)
517: dof[i*dmmoab->numFields+field] = dmmoab->lidmap[(PetscInt)points[i]]*dmmoab->numFields+field;
518: }
519: }
520: else {
521: /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
522: /* TODO: eliminate the limitation using PetscSection to manage DOFs */
523: for (field=0; field<dmmoab->numFields; ++field) {
524: offset = field*dmmoab->n;
525: for (i=0; i<npoints; ++i)
526: dof[i*dmmoab->numFields+field] = dmmoab->lidmap[(PetscInt)points[i]]+offset;
527: }
528: }
529: return(0);
530: }
535: /*@C
536: DMMoabGetDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
537: array of MOAB EntityHandles. It is useful when performing 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 global degree-of-freedom index array in the discrete representation (Vec, Mat)
550: Level: intermediate
552: .keywords: discretization manager, get, global degrees of freedom
553: 554: .seealso: DMMoabGetDofsLocal(), DMMoabGetDofs(), DMMoabGetDofsBlockedLocal()
555: @*/
556: PetscErrorCodeDMMoabGetDofsBlocked(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->gidmap[(PetscInt)points[i]];
573: }
574: return(0);
575: }
580: /*@C
581: DMMoabGetDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
582: array of MOAB EntityHandles. It is useful when performing local Blocked(Get/Set) methods in computation
583: of element residuals and assembly of the discrete systems when all fields are co-located.
585: Not Collective
587: Input Parameters:
588: + dm - the discretization manager object
589: . npoints - the total number of Entities in the points array
590: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
592: Output Parameter:
593: + dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat)
595: Level: intermediate
597: .keywords: discretization manager, get, global degrees of freedom
598: 599: .seealso: DMMoabGetDofsLocal(), DMMoabGetDofs(), DMMoabGetDofsBlockedLocal()
600: @*/
601: PetscErrorCodeDMMoabGetDofsBlockedLocal(DM dm,PetscInt npoints,const moab::EntityHandle* points,PetscInt* dof)
602: {
603: PetscInt i;
604: DM_Moab *dmmoab;
605: PetscErrorCode ierr;
610: dmmoab = (DM_Moab*)(dm)->data;
612: if (!dof) {
613: PetscMalloc1(npoints,&dof);
614: }
616: for (i=0; i<npoints; ++i)
617: dof[i] = dmmoab->lidmap[(PetscInt)points[i]];
618: return(0);
619: }
624: /*@C
625: DMMoabGetVertexDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
626: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
627: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
629: Not Collective
631: Input Parameters:
632: + dm - the discretization manager object
634: Output Parameter:
635: + dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
637: Level: intermediate
639: .keywords: discretization manager, get, blocked degrees of freedom
640: 641: .seealso: DMMoabGetVertexDofsBlockedLocal(), DMMoabGetDofsBlocked(), DMMoabGetDofsBlockedLocal()
642: @*/
643: PetscErrorCodeDMMoabGetVertexDofsBlocked(DM dm,PetscInt** dof)644: {
645: DM_Moab *dmmoab;
649: dmmoab = (DM_Moab*)(dm)->data;
651: *dof = dmmoab->gidmap;
652: return(0);
653: }
658: /*@C
659: DMMoabGetVertexDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
660: array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
661: where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.
663: Not Collective
665: Input Parameters:
666: + dm - the discretization manager object
668: Output Parameter:
669: + dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering
671: Level: intermediate
673: .keywords: discretization manager, get, blocked degrees of freedom
674: 675: .seealso: DMMoabGetVertexDofsBlocked(), DMMoabGetDofsBlocked(), DMMoabGetDofsBlockedLocal()
676: @*/
677: PetscErrorCodeDMMoabGetVertexDofsBlockedLocal(DM dm,PetscInt** dof)678: {
679: DM_Moab *dmmoab;
684: dmmoab = (DM_Moab*)(dm)->data;
686: *dof = dmmoab->lidmap;
687: return(0);
688: }