Actual source code: dmmbfield.cxx

petsc-3.14.6 2021-03-30
Report Typos and Errors
  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: }