Actual source code: dmfield.c

  1: #include <petsc/private/dmfieldimpl.h>
  2: #include <petsc/private/petscfeimpl.h>
  3: #include <petscdmplex.h>

  5: const char *const DMFieldContinuities[] = {"VERTEX", "EDGE", "FACET", "CELL", NULL};

  7: PETSC_INTERN PetscErrorCode DMFieldCreate(DM dm, PetscInt numComponents, DMFieldContinuity continuity, DMField *field)
  8: {
  9:   DMField b;

 13:   DMFieldInitializePackage();

 15:   PetscHeaderCreate(b, DMFIELD_CLASSID, "DMField", "Field over DM", "DM", PetscObjectComm((PetscObject)dm), DMFieldDestroy, DMFieldView);
 16:   PetscObjectReference((PetscObject)dm);
 17:   b->dm            = dm;
 18:   b->continuity    = continuity;
 19:   b->numComponents = numComponents;
 20:   *field           = b;
 21:   return 0;
 22: }

 24: /*@
 25:    DMFieldDestroy - destroy a `DMField`

 27:    Collective

 29:    Input Parameter:
 30: .  field - address of `DMField`

 32:    Level: advanced

 34: .seealso: `DMField`, `DMFieldCreate()`
 35: @*/
 36: PetscErrorCode DMFieldDestroy(DMField *field)
 37: {
 38:   if (!*field) return 0;
 40:   if (--((PetscObject)(*field))->refct > 0) {
 41:     *field = NULL;
 42:     return 0;
 43:   }
 44:   PetscTryTypeMethod((*field), destroy);
 45:   DMDestroy(&((*field)->dm));
 46:   PetscHeaderDestroy(field);
 47:   return 0;
 48: }

 50: /*@C
 51:    DMFieldView - view a DMField

 53:    Collective

 55:    Input Parameters:
 56: +  field - DMField
 57: -  viewer - viewer to display field, for example PETSC_VIEWER_STDOUT_WORLD

 59:    Level: advanced

 61: .seealso: `DMField`, `DMFieldCreate()`
 62: @*/
 63: PetscErrorCode DMFieldView(DMField field, PetscViewer viewer)
 64: {
 65:   PetscBool iascii;

 68:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)field), &viewer);
 71:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
 72:   if (iascii) {
 73:     PetscObjectPrintClassNamePrefixType((PetscObject)field, viewer);
 74:     PetscViewerASCIIPushTab(viewer);
 75:     PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " components\n", field->numComponents);
 76:     PetscViewerASCIIPrintf(viewer, "%s continuity\n", DMFieldContinuities[field->continuity]);
 77:     PetscViewerPushFormat(viewer, PETSC_VIEWER_DEFAULT);
 78:     DMView(field->dm, viewer);
 79:     PetscViewerPopFormat(viewer);
 80:   }
 81:   PetscTryTypeMethod(field, view, viewer);
 82:   if (iascii) PetscViewerASCIIPopTab(viewer);
 83:   return 0;
 84: }

 86: /*@C
 87:    DMFieldSetType - set the `DMField` implementation

 89:    Collective on field

 91:    Input Parameters:
 92: +  field - the `DMField` context
 93: -  type - a known method

 95:    Notes:
 96:    See "include/petscvec.h" for available methods (for instance)
 97: +    `DMFIELDDA`    - a field defined only by its values at the corners of a `DMDA`
 98: .    `DMFIELDDS`    - a field defined by a discretization over a mesh set with `DMSetField()`
 99: -    `DMFIELDSHELL` - a field defined by arbitrary callbacks

101:   Level: advanced

103: .seealso: `DMField`, `DMFieldGetType()`, `DMFieldType`,
104: @*/
105: PetscErrorCode DMFieldSetType(DMField field, DMFieldType type)
106: {
107:   PetscBool match;
108:   PetscErrorCode (*r)(DMField);


113:   PetscObjectTypeCompare((PetscObject)field, type, &match);
114:   if (match) return 0;

116:   PetscFunctionListFind(DMFieldList, type, &r);
118:   /* Destroy the previous private DMField context */
119:   PetscTryTypeMethod(field, destroy);

121:   PetscMemzero(field->ops, sizeof(*field->ops));
122:   PetscObjectChangeTypeName((PetscObject)field, type);
123:   field->ops->create = r;
124:   (*r)(field);
125:   return 0;
126: }

128: /*@C
129:   DMFieldGetType - Gets the `DMFieldType` name (as a string) from the `DMField`.

131:   Not Collective

133:   Input Parameter:
134: . field  - The `DMField` context

136:   Output Parameter:
137: . type - The `DMFieldType` name

139:   Level: advanced

141: .seealso: `DMField`, `DMFieldSetType()`, `DMFieldType`
142: @*/
143: PetscErrorCode DMFieldGetType(DMField field, DMFieldType *type)
144: {
147:   DMFieldRegisterAll();
148:   *type = ((PetscObject)field)->type_name;
149:   return 0;
150: }

152: /*@
153:   DMFieldGetNumComponents - Returns the number of components in the field

155:   Not collective

157:   Input Parameter:
158: . field - The `DMField` object

160:   Output Parameter:
161: . nc - The number of field components

163:   Level: intermediate

165: .seealso: `DMField`, `DMFieldEvaluate()`
166: @*/
167: PetscErrorCode DMFieldGetNumComponents(DMField field, PetscInt *nc)
168: {
171:   *nc = field->numComponents;
172:   return 0;
173: }

175: /*@
176:   DMFieldGetDM - Returns the `DM` for the manifold over which the field is defined.

178:   Not collective

180:   Input Parameter:
181: . field - The `DMField` object

183:   Output Parameter:
184: . dm - The `DM` object

186:   Level: intermediate

188: .seealso: `DMField`, `DM`, `DMFieldEvaluate()`
189: @*/
190: PetscErrorCode DMFieldGetDM(DMField field, DM *dm)
191: {
194:   *dm = field->dm;
195:   return 0;
196: }

198: /*@
199:   DMFieldEvaluate - Evaluate the field and its derivatives on a set of points

201:   Collective on points

203:   Input Parameters:
204: + field - The `DMField` object
205: . points - The points at which to evaluate the field.  Should have size d x n,
206:            where d is the coordinate dimension of the manifold and n is the number
207:            of points
208: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
209:              If the field is complex and datatype is `PETSC_REAL`, the real part of the
210:              field is returned.

212:   Output Parameters:
213: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
214:       If B is not NULL, the values of the field are written in this array, varying first by component,
215:       then by point.
216: . D - pointer to data of size d * c * n * sizeof(datatype).
217:       If D is not NULL, the values of the field's spatial derivatives are written in this array,
218:       varying first by the partial derivative component, then by field component, then by point.
219: - H - pointer to data of size d * d * c * n * sizeof(datatype).
220:       If H is not NULL, the values of the field's second spatial derivatives are written in this array,
221:       varying first by the second partial derivative component, then by field component, then by point.

223:   Level: intermediate

225: .seealso: `DMField`, `DMFieldGetDM()`, `DMFieldGetNumComponents()`, `DMFieldEvaluateFE()`, `DMFieldEvaluateFV()`, `PetscDataType`
226: @*/
227: PetscErrorCode DMFieldEvaluate(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
228: {
234:   if (field->ops->evaluate) {
235:     (*field->ops->evaluate)(field, points, datatype, B, D, H);
236:   } else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented for this type");
237:   return 0;
238: }

240: /*@
241:   DMFieldEvaluateFE - Evaluate the field and its derivatives on a set of points mapped from
242:   quadrature points on a reference point.  The derivatives are taken with respect to the
243:   reference coordinates.

245:   Not collective

247:   Input Parameters:
248: + field - The `DMField` object
249: . cellIS - Index set for cells on which to evaluate the field
250: . points - The quadature containing the points in the reference cell at which to evaluate the field.
251: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
252:              If the field is complex and datatype is PETSC_REAL, the real part of the
253:              field is returned.

255:   Output Parameters:
256: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
257:       If B is not NULL, the values of the field are written in this array, varying first by component,
258:       then by point.
259: . D - pointer to data of size d * c * n * sizeof(datatype).
260:       If D is not NULL, the values of the field's spatial derivatives are written in this array,
261:       varying first by the partial derivative component, then by field component, then by point.
262: - H - pointer to data of size d * d * c * n * sizeof(datatype).
263:       If H is not NULL, the values of the field's second spatial derivatives are written in this array,
264:       varying first by the second partial derivative component, then by field component, then by point.

266:   Level: intermediate

268: .seealso: `DMField`, `DM`, `DMFieldGetNumComponents()`, `DMFieldEvaluate()`, `DMFieldEvaluateFV()`
269: @*/
270: PetscErrorCode DMFieldEvaluateFE(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
271: {
278:   if (field->ops->evaluateFE) {
279:     (*field->ops->evaluateFE)(field, cellIS, points, datatype, B, D, H);
280:   } else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented for this type");
281:   return 0;
282: }

284: /*@
285:   DMFieldEvaluateFV - Evaluate the mean of a field and its finite volume derivatives on a set of points.

287:   Not collective

289:   Input Parameters:
290: + field - The `DMField` object
291: . cellIS - Index set for cells on which to evaluate the field
292: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
293:              If the field is complex and datatype is `PETSC_REAL`, the real part of the
294:              field is returned.

296:   Output Parameters:
297: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
298:       If B is not NULL, the values of the field are written in this array, varying first by component,
299:       then by point.
300: . D - pointer to data of size d * c * n * sizeof(datatype).
301:       If D is not NULL, the values of the field's spatial derivatives are written in this array,
302:       varying first by the partial derivative component, then by field component, then by point.
303: - H - pointer to data of size d * d * c * n * sizeof(datatype).
304:       If H is not NULL, the values of the field's second spatial derivatives are written in this array,
305:       varying first by the second partial derivative component, then by field component, then by point.

307:   Level: intermediate

309: .seealso: `DMField`, `IS`, `DMFieldGetNumComponents()`, `DMFieldEvaluate()`, `DMFieldEvaluateFE()`, `PetscDataType`
310: @*/
311: PetscErrorCode DMFieldEvaluateFV(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
312: {
318:   if (field->ops->evaluateFV) {
319:     (*field->ops->evaluateFV)(field, cellIS, datatype, B, D, H);
320:   } else SETERRQ(PetscObjectComm((PetscObject)field), PETSC_ERR_SUP, "Not implemented for this type");
321:   return 0;
322: }

324: /*@
325:   DMFieldGetDegree - Get the polynomial degree of a field when pulled back onto the
326:   reference element

328:   Not collective

330:   Input Parameters:
331: + field - the `DMField` object
332: - cellIS - the index set of points over which we want know the invariance

334:   Output Parameters:
335: + minDegree - the degree of the largest polynomial space contained in the field on each element
336: - maxDegree - the largest degree of the smallest polynomial space containing the field on any element

338:   Level: intermediate

340: .seealso: `DMField`, `IS`, `DMFieldEvaluateFE()`
341: @*/
342: PetscErrorCode DMFieldGetDegree(DMField field, IS cellIS, PetscInt *minDegree, PetscInt *maxDegree)
343: {

349:   if (minDegree) *minDegree = -1;
350:   if (maxDegree) *maxDegree = PETSC_MAX_INT;

352:   if (field->ops->getDegree) (*field->ops->getDegree)(field, cellIS, minDegree, maxDegree);
353:   return 0;
354: }

356: /*@
357:   DMFieldCreateDefaultQuadrature - Creates a quadrature sufficient to integrate the field on the selected
358:   points via pullback onto the reference element

360:   Not collective

362:   Input Parameters:
363: + field - the `DMField` object
364: - pointIS - the index set of points over which we wish to integrate the field

366:   Output Parameter:
367: . quad - a `PetscQuadrature` object

369:   Level: developer

371: .seealso: `DMField`, `PetscQuadrature`, `IS`, `DMFieldEvaluteFE()`, `DMFieldGetDegree()`
372: @*/
373: PetscErrorCode DMFieldCreateDefaultQuadrature(DMField field, IS pointIS, PetscQuadrature *quad)
374: {

379:   *quad = NULL;
380:   PetscTryTypeMethod(field, createDefaultQuadrature, pointIS, quad);
381:   return 0;
382: }

384: /*@C
385:   DMFieldCreateFEGeom - Compute and create the geometric factors of a coordinate field

387:   Not collective

389:   Input Parameters:
390: + field - the `DMField` object
391: . pointIS - the index set of points over which we wish to integrate the field
392: . quad - the quadrature points at which to evaluate the geometric factors
393: - faceData - whether additional data for facets (the normal vectors and adjacent cells) should
394:   be calculated

396:   Output Parameter:
397: . geom - the geometric factors

399:   Level: developer

401: .seealso: `DMField`, `PetscQuadrature`, `IS`, `PetscFEGeom`, `DMFieldEvaluateFE()`, `DMFieldCreateDefaulteQuadrature()`, `DMFieldGetDegree()`
402: @*/
403: PetscErrorCode DMFieldCreateFEGeom(DMField field, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
404: {
405:   PetscInt     dim, dE;
406:   PetscInt     nPoints;
407:   PetscInt     maxDegree;
408:   PetscFEGeom *g;

413:   ISGetLocalSize(pointIS, &nPoints);
414:   dE = field->numComponents;
415:   PetscFEGeomCreate(quad, nPoints, dE, faceData, &g);
416:   DMFieldEvaluateFE(field, pointIS, quad, PETSC_REAL, g->v, g->J, NULL);
417:   dim = g->dim;
418:   if (dE > dim) {
419:     /* space out J and make square Jacobians */
420:     PetscInt i, j, k, N = g->numPoints * g->numCells;

422:     for (i = N - 1; i >= 0; i--) {
423:       PetscReal J[9];

425:       for (j = 0; j < dE; j++) {
426:         for (k = 0; k < dim; k++) J[j * dE + k] = g->J[i * dE * dim + j * dim + k];
427:       }
428:       switch (dim) {
429:       case 0:
430:         for (j = 0; j < dE; j++) {
431:           for (k = 0; k < dE; k++) J[j * dE + k] = (j == k) ? 1. : 0.;
432:         }
433:         break;
434:       case 1:
435:         if (dE == 2) {
436:           PetscReal norm = PetscSqrtReal(J[0] * J[0] + J[2] * J[2]);

438:           J[1] = -J[2] / norm;
439:           J[3] = J[0] / norm;
440:         } else {
441:           PetscReal inorm = 1. / PetscSqrtReal(J[0] * J[0] + J[3] * J[3] + J[6] * J[6]);
442:           PetscReal x     = J[0] * inorm;
443:           PetscReal y     = J[3] * inorm;
444:           PetscReal z     = J[6] * inorm;

446:           if (x > 0.) {
447:             PetscReal inv1pX = 1. / (1. + x);

449:             J[1] = -y;
450:             J[2] = -z;
451:             J[4] = 1. - y * y * inv1pX;
452:             J[5] = -y * z * inv1pX;
453:             J[7] = -y * z * inv1pX;
454:             J[8] = 1. - z * z * inv1pX;
455:           } else {
456:             PetscReal inv1mX = 1. / (1. - x);

458:             J[1] = z;
459:             J[2] = y;
460:             J[4] = -y * z * inv1mX;
461:             J[5] = 1. - y * y * inv1mX;
462:             J[7] = 1. - z * z * inv1mX;
463:             J[8] = -y * z * inv1mX;
464:           }
465:         }
466:         break;
467:       case 2: {
468:         PetscReal inorm;

470:         J[2] = J[3] * J[7] - J[6] * J[4];
471:         J[5] = J[6] * J[1] - J[0] * J[7];
472:         J[8] = J[0] * J[4] - J[3] * J[1];

474:         inorm = 1. / PetscSqrtReal(J[2] * J[2] + J[5] * J[5] + J[8] * J[8]);

476:         J[2] *= inorm;
477:         J[5] *= inorm;
478:         J[8] *= inorm;
479:       } break;
480:       }
481:       for (j = 0; j < dE * dE; j++) g->J[i * dE * dE + j] = J[j];
482:     }
483:   }
484:   PetscFEGeomComplete(g);
485:   DMFieldGetDegree(field, pointIS, NULL, &maxDegree);
486:   g->isAffine = (maxDegree <= 1) ? PETSC_TRUE : PETSC_FALSE;
487:   if (faceData) (*field->ops->computeFaceData)(field, pointIS, quad, g);
488:   *geom = g;
489:   return 0;
490: }