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:   std::string        tag_name;

 30:   dmmoab = (DM_Moab *)(dm)->data;


 34:   /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
 35:   merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
 36:   MBERRNM(merr);

 38:   DMMoabGetVecTag(fvec, &vtag);

 40:   merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
 41:   if (!tag_name.length() && merr != moab::MB_SUCCESS) {
 42:     VecGetArrayRead(fvec, &varray);
 43:     /* use the entity handle and the Dof index to set the right value */
 44:     merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)varray);
 45:     MBERRNM(merr);
 46:     VecRestoreArrayRead(fvec, &varray);
 47:   } else {
 48:     PetscMalloc1(dmmoab->nloc, &farray);
 49:     /* we are using a MOAB Vec - directly copy the tag data to new one */
 50:     merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void *)farray);
 51:     MBERRNM(merr);
 52:     merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)farray);
 53:     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);
 59:   MBERRNM(merr);
 60: #endif
 61:   return 0;
 62: }

 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:   PetscInt              i, ifield;
 88:   std::string           tag_name;
 89:   moab::Range::iterator iter;

 92:   dmmoab = (DM_Moab *)(dm)->data;

 94:   /* get the Tag corresponding to the global vector - possible that there is no tag associated.. */
 95:   DMMoabGetVecTag(fvec, &vtag);
 96:   merr = dmmoab->mbiface->tag_get_name(vtag, tag_name);
 97:   PetscMalloc1(dmmoab->nloc, &farray);
 98:   if (!tag_name.length() && merr != moab::MB_SUCCESS) {
 99:     /* not a MOAB vector - use VecGetSubVector to get the parts as needed */
100:     VecGetArrayRead(fvec, &rarray);
101:     for (ifield = 0; ifield < dmmoab->numFields; ++ifield) {
102:       /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
103:       merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
104:       MBERRNM(merr);

106:       for (i = 0; i < dmmoab->nloc; i++) farray[i] = (dmmoab->bs == 1 ? rarray[ifield * dmmoab->nloc + i] : rarray[i * dmmoab->numFields + ifield]);

108:       /* use the entity handle and the Dof index to set the right value */
109:       merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)farray);
110:       MBERRNM(merr);
111:     }
112:     VecRestoreArrayRead(fvec, &rarray);
113:   } else {
114:     PetscMalloc1(dmmoab->nloc * dmmoab->numFields, &varray);

116:     /* we are using a MOAB Vec - directly copy the tag data to new one */
117:     merr = dmmoab->mbiface->tag_get_data(vtag, *dmmoab->vowned, (void *)varray);
118:     MBERRNM(merr);
119:     for (ifield = 0; ifield < dmmoab->numFields; ++ifield) {
120:       /* Create a tag in MOAB mesh to index and keep track of number of Petsc vec tags */
121:       merr = dmmoab->mbiface->tag_get_handle(dmmoab->fieldNames[ifield], 1, moab::MB_TYPE_DOUBLE, ntag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT);
122:       MBERRNM(merr);

124:       /* we are using a MOAB Vec - directly copy the tag data to new one */
125:       for (i = 0; i < dmmoab->nloc; i++) farray[i] = (dmmoab->bs == 1 ? varray[ifield * dmmoab->nloc + i] : varray[i * dmmoab->numFields + ifield]);

127:       merr = dmmoab->mbiface->tag_set_data(ntag, *dmmoab->vowned, (const void *)farray);
128:       MBERRNM(merr);

130: #ifdef MOAB_HAVE_MPI
131:       /* make sure the parallel exchange for ghosts are done appropriately */
132:       merr = dmmoab->pcomm->exchange_tags(ntag, *dmmoab->vlocal);
133:       MBERRNM(merr);
134: #endif
135:     }
136:     PetscFree(varray);
137:   }
138:   PetscFree(farray);
139:   return 0;
140: }

142: /*@C
143:   DMMoabSetFieldNames - Sets the number of fields and their names to be managed by the DM

145:   Not Collective

147:   Input Parameters:
148: + dm     - the discretization manager object
149: . numFields - the total number of fields
150: - fields - the array containing the names of each field (component); Can be NULL.

152:   Level: intermediate

154: .seealso: `DMMoabGetFieldName()`, `DMMoabSetFieldName()`
155: @*/
156: PetscErrorCode DMMoabSetFieldNames(DM dm, PetscInt numFields, const char *fields[])
157: {
158:   PetscInt i;
159:   DM_Moab *dmmoab;

162:   dmmoab = (DM_Moab *)(dm)->data;

164:   /* first deallocate the existing field structure */
165:   if (dmmoab->fieldNames) {
166:     for (i = 0; i < dmmoab->numFields; i++) PetscFree(dmmoab->fieldNames[i]);
167:     PetscFree(dmmoab->fieldNames);
168:   }

170:   /* now re-allocate and assign field names  */
171:   dmmoab->numFields = numFields;
172:   PetscMalloc1(numFields, &dmmoab->fieldNames);
173:   if (fields) {
174:     for (i = 0; i < dmmoab->numFields; i++) PetscStrallocpy(fields[i], (char **)&dmmoab->fieldNames[i]);
175:   }
176:   DMSetNumFields(dm, numFields);
177:   return 0;
178: }

180: /*@C
181:   DMMoabGetFieldName - Gets the names of individual field components in multicomponent
182:   vectors associated with a DMDA.

184:   Not Collective

186:   Input Parameters:
187: + dm     - the discretization manager object
188: - field - field number for the DMMoab (0, 1, ... dof-1), where dof indicates the
189:         number of degrees of freedom per node within the DMMoab

191:   Output Parameter:
192: . fieldName - the name of the field (component)

194:   Level: intermediate

196: .seealso: `DMMoabSetFieldName()`, `DMMoabSetFields()`
197: @*/
198: PetscErrorCode DMMoabGetFieldName(DM dm, PetscInt field, const char **fieldName)
199: {
200:   DM_Moab *dmmoab;

203:   dmmoab = (DM_Moab *)(dm)->data;

206:   *fieldName = dmmoab->fieldNames[field];
207:   return 0;
208: }

210: /*@C
211:   DMMoabSetFieldName - Sets the name of a field (component) managed by the DM

213:   Not Collective

215:   Input Parameters:
216: + dm     - the discretization manager object
217: . field - the field number
218: - fieldName - the field (component) name

220:   Level: intermediate
221:   Notes:
222:     Can only be called after DMMoabSetFields supplied with correct numFields

224: .seealso: `DMMoabGetFieldName()`, `DMMoabSetFields()`
225: @*/
226: PetscErrorCode DMMoabSetFieldName(DM dm, PetscInt field, const char *fieldName)
227: {
228:   DM_Moab *dmmoab;


233:   dmmoab = (DM_Moab *)(dm)->data;

236:   if (dmmoab->fieldNames[field]) PetscFree(dmmoab->fieldNames[field]);
237:   PetscStrallocpy(fieldName, (char **)&dmmoab->fieldNames[field]);
238:   return 0;
239: }

241: /*@C
242:   DMMoabGetFieldDof - Gets the global degree-of-freedom of a field (component) defined on a
243:   particular MOAB EntityHandle.

245:   Not Collective

247:   Input Parameters:
248: + dm     - the discretization manager object
249: . point - the MOAB EntityHandle container which holds the field degree-of-freedom values
250: - field - the field (component) index

252:   Output Parameter:
253: . dof - the global degree-of-freedom index corresponding to the field in the discrete representation (Vec, Mat)

255:   Level: beginner

257: .seealso: `DMMoabGetFieldDofs()`, `DMMoabGetFieldDofsLocal()`
258: @*/
259: PetscErrorCode DMMoabGetFieldDof(DM dm, moab::EntityHandle point, PetscInt field, PetscInt *dof)
260: {
261:   DM_Moab *dmmoab;

264:   dmmoab = (DM_Moab *)(dm)->data;

266:   *dof = (dmmoab->bs == 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(point) - dmmoab->seqstart] + field * dmmoab->n : dmmoab->gidmap[dmmoab->mbiface->id_from_handle(point) - dmmoab->seqstart] * dmmoab->numFields + field);
267:   return 0;
268: }

270: /*@C
271:   DMMoabGetFieldDofs - Gets the global degree-of-freedom of a field (component) defined on an
272:   array of MOAB EntityHandles.

274:   Not Collective

276:   Input Parameters:
277: + dm     - the discretization manager object
278: . npoints - the total number of Entities in the points array
279: . points - the MOAB EntityHandle container array 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 array corresponding to the field in the discrete representation (Vec, Mat)

285:   Level: intermediate

287: .seealso: `DMMoabGetFieldDof()`, `DMMoabGetFieldDofsLocal()`
288: @*/
289: PetscErrorCode DMMoabGetFieldDofs(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt field, PetscInt *dof)
290: {
291:   PetscInt i;
292:   DM_Moab *dmmoab;

296:   dmmoab = (DM_Moab *)(dm)->data;

298:   if (!dof) PetscMalloc1(npoints, &dof);

300:   /* compute the DOF based on local blocking in the fields */
301:   /* We also assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
302:   /* TODO: eliminate the limitation using PetscSection to manage DOFs */
303:   for (i = 0; i < npoints; ++i)
304:     dof[i] = (dmmoab->bs == 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + field * dmmoab->n : dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field);
305:   return 0;
306: }

308: /*@C
309:   DMMoabGetFieldDofsLocal - Gets the local degrees-of-freedom of a field (component) defined on an
310:   array of MOAB EntityHandles.

312:   Not Collective

314:   Input Parameters:
315: + dm     - the discretization manager object
316: . npoints - the total number of Entities in the points array
317: . points - the MOAB EntityHandle container array which holds the field degree-of-freedom values
318: - field - the field (component) index

320:   Output Parameter:
321: . dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)

323:   Level: intermediate

325: .seealso: `DMMoabGetFieldDof()`, `DMMoabGetFieldDofs()`
326: @*/
327: PetscErrorCode DMMoabGetFieldDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt field, PetscInt *dof)
328: {
329:   PetscInt i;
330:   DM_Moab *dmmoab;

334:   dmmoab = (DM_Moab *)(dm)->data;

336:   if (!dof) PetscMalloc1(npoints, &dof);

338:   /* compute the DOF based on local blocking in the fields */
339:   /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
340:   /* TODO: eliminate the limitation using PetscSection to manage DOFs */
341:   for (i = 0; i < npoints; ++i) {
342:     dof[i] = (dmmoab->bs > 1 ? dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field : dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + field * dmmoab->n);
343:   }
344:   return 0;
345: }

347: /*@C
348:   DMMoabGetDofs - Gets the global degree-of-freedom for all fields (components) defined on an
349:   array of MOAB EntityHandles.

351:   Not Collective

353:   Input Parameters:
354: + dm     - the discretization manager object
355: . npoints - the total number of Entities in the points array
356: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values

358:   Output Parameter:
359: . dof - the global degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)

361:   Level: intermediate

363: .seealso: `DMMoabGetFieldDofs()`, `DMMoabGetDofsLocal()`, `DMMoabGetDofsBlocked()`
364: @*/
365: PetscErrorCode DMMoabGetDofs(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
366: {
367:   PetscInt i, field, offset;
368:   DM_Moab *dmmoab;

372:   dmmoab = (DM_Moab *)(dm)->data;

374:   if (!dof) PetscMalloc1(dmmoab->numFields * npoints, &dof);

376:   /* compute the DOF based on local blocking in the fields */
377:   /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
378:   /* TODO: eliminate the limitation using PetscSection to manage DOFs */
379:   for (field = 0; field < dmmoab->numFields; ++field) {
380:     offset = field * dmmoab->n;
381:     for (i = 0; i < npoints; ++i)
382:       dof[i * dmmoab->numFields + field] = (dmmoab->bs > 1 ? dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field : dmmoab->gidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + offset);
383:   }
384:   return 0;
385: }

387: /*@C
388:   DMMoabGetDofsLocal - Gets the local degree-of-freedom for all fields (components) defined on an
389:   array of MOAB EntityHandles.

391:   Not Collective

393:   Input Parameters:
394: + dm     - the discretization manager object
395: . npoints - the total number of Entities in the points array
396: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values

398:   Output Parameter:
399: . dof - the local degree-of-freedom index array corresponding to the field in the discrete representation (Vec, Mat)

401:   Level: intermediate

403: .seealso: `DMMoabGetFieldDofs()`, `DMMoabGetDofs()`, `DMMoabGetDofsBlocked()`
404: @*/
405: PetscErrorCode DMMoabGetDofsLocal(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
406: {
407:   PetscInt i, field, offset;
408:   DM_Moab *dmmoab;

412:   dmmoab = (DM_Moab *)(dm)->data;

414:   if (!dof) PetscMalloc1(dmmoab->numFields * npoints, &dof);

416:   /* compute the DOF based on local blocking in the fields */
417:   /* assume all fields have equal distribution; i.e., all fields are either defined on vertices or elements and not on a mixture */
418:   /* TODO: eliminate the limitation using PetscSection to manage DOFs */
419:   for (field = 0; field < dmmoab->numFields; ++field) {
420:     offset = field * dmmoab->n;
421:     for (i = 0; i < npoints; ++i)
422:       dof[i * dmmoab->numFields + field] = (dmmoab->bs > 1 ? dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] * dmmoab->numFields + field : dmmoab->lidmap[dmmoab->mbiface->id_from_handle(points[i]) - dmmoab->seqstart] + offset);
423:   }
424:   return 0;
425: }

427: /*@C
428:   DMMoabGetDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
429:   array of MOAB EntityHandles. It is useful when performing Blocked(Get/Set) methods in computation
430:   of element residuals and assembly of the discrete systems when all fields are co-located.

432:   Not Collective

434:   Input Parameters:
435: + dm     - the discretization manager object
436: . npoints - the total number of Entities in the points array
437: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values

439:   Output Parameter:
440: . dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat)

442:   Level: intermediate

444: .seealso: `DMMoabGetDofsLocal()`, `DMMoabGetDofs()`, `DMMoabGetDofsBlockedLocal()`
445: @*/
446: PetscErrorCode DMMoabGetDofsBlocked(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
447: {
448:   PetscInt i;
449:   DM_Moab *dmmoab;

453:   dmmoab = (DM_Moab *)(dm)->data;

455:   if (!dof) PetscMalloc1(npoints, &dof);

457:   for (i = 0; i < npoints; ++i) dof[i] = dmmoab->gidmap[(PetscInt)points[i] - dmmoab->seqstart];
458:   return 0;
459: }

461: /*@C
462:   DMMoabGetDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
463:   array of MOAB EntityHandles. It is useful when performing local Blocked(Get/Set) methods in computation
464:   of element residuals and assembly of the discrete systems when all fields are co-located.

466:   Not Collective

468:   Input Parameters:
469: + dm     - the discretization manager object
470: . npoints - the total number of Entities in the points array
471: - points - the MOAB EntityHandle container array which holds the field degree-of-freedom values

473:   Output Parameter:
474: . dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat)

476:   Level: intermediate

478: .seealso: `DMMoabGetDofsLocal()`, `DMMoabGetDofs()`, `DMMoabGetDofsBlockedLocal()`
479: @*/
480: PetscErrorCode DMMoabGetDofsBlockedLocal(DM dm, PetscInt npoints, const moab::EntityHandle *points, PetscInt *dof)
481: {
482:   PetscInt i;
483:   DM_Moab *dmmoab;

487:   dmmoab = (DM_Moab *)(dm)->data;

489:   if (!dof) PetscMalloc1(npoints, &dof);

491:   for (i = 0; i < npoints; ++i) dof[i] = dmmoab->lidmap[(PetscInt)points[i] - dmmoab->seqstart];
492:   return 0;
493: }

495: /*@C
496:   DMMoabGetVertexDofsBlocked - Gets the global degree-of-freedom for the first field (component) defined on an
497:   array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
498:   where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.

500:   Not Collective

502:   Input Parameters:
503: . dm     - the discretization manager object

505:   Output Parameter:
506: . dof - the blocked global degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering

508:   Level: intermediate

510: .seealso: `DMMoabGetVertexDofsBlockedLocal()`, `DMMoabGetDofsBlocked()`, `DMMoabGetDofsBlockedLocal()`
511: @*/
512: PetscErrorCode DMMoabGetVertexDofsBlocked(DM dm, PetscInt **dof)
513: {
514:   DM_Moab *dmmoab;

517:   dmmoab = (DM_Moab *)(dm)->data;

519:   *dof = dmmoab->gidmap;
520:   return 0;
521: }

523: /*@C
524:   DMMoabGetVertexDofsBlockedLocal - Gets the local degree-of-freedom for the first field (component) defined on an
525:   array of locally owned MOAB mesh vertices. It's utility is when performing Finite-Difference type calculations
526:   where vertex traversal is faster than element-wise assembly that is typically done in FEM calculations.

528:   Not Collective

530:   Input Parameters:
531: . dm     - the discretization manager object

533:   Output Parameter:
534: . dof - the blocked local degree-of-freedom index array in the discrete representation (Vec, Mat) that is vertex-based based on local numbering

536:   Level: intermediate

538: .seealso: `DMMoabGetVertexDofsBlocked()`, `DMMoabGetDofsBlocked()`, `DMMoabGetDofsBlockedLocal()`
539: @*/
540: PetscErrorCode DMMoabGetVertexDofsBlockedLocal(DM dm, PetscInt **dof)
541: {
542:   DM_Moab *dmmoab;

546:   dmmoab = (DM_Moab *)(dm)->data;

548:   *dof = dmmoab->lidmap;
549:   return 0;
550: }