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;

     11:   PetscFunctionBegin;
     13:   PetscAssertPointer(field, 4);
     14:   PetscCall(DMFieldInitializePackage());

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

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

     28:   Collective

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

     33:   Level: advanced

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

     52: /*@
     53:   DMFieldView - view a `DMField`

     55:   Collective

     57:   Input Parameters:
     58: + field  - `DMField`
     59: - viewer - viewer to display field, for example `PETSC_VIEWER_STDOUT_WORLD`

     61:   Level: advanced

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

     69:   PetscFunctionBegin;
     71:   if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)field), &viewer));
     73:   PetscCheckSameComm(field, 1, viewer, 2);
     74:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
     75:   if (iascii) {
     76:     PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)field, viewer));
     77:     PetscCall(PetscViewerASCIIPushTab(viewer));
     78:     PetscCall(PetscViewerASCIIPrintf(viewer, "%" PetscInt_FMT " components\n", field->numComponents));
     79:     PetscCall(PetscViewerASCIIPrintf(viewer, "%s continuity\n", DMFieldContinuities[field->continuity]));
     80:     PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_DEFAULT));
     81:     PetscCall(DMView(field->dm, viewer));
     82:     PetscCall(PetscViewerPopFormat(viewer));
     83:   }
     84:   PetscTryTypeMethod(field, view, viewer);
     85:   if (iascii) PetscCall(PetscViewerASCIIPopTab(viewer));
     86:   PetscFunctionReturn(PETSC_SUCCESS);
     87: }

     89: /*@
     90:   DMFieldSetType - set the `DMField` implementation

     92:   Collective

     94:   Input Parameters:
     95: + field - the `DMField` context
     96: - type  - a known method

     98:   Level: advanced

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

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

    113:   PetscFunctionBegin;
    115:   PetscAssertPointer(type, 2);

    117:   PetscCall(PetscObjectTypeCompare((PetscObject)field, type, &match));
    118:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

    120:   PetscCall(PetscFunctionListFind(DMFieldList, type, &r));
    121:   PetscCheck(r, PetscObjectComm((PetscObject)field), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unable to find requested DMField type %s", type);
    122:   /* Destroy the previous private DMField context */
    123:   PetscTryTypeMethod(field, destroy);

    125:   PetscCall(PetscMemzero(field->ops, sizeof(*field->ops)));
    126:   PetscCall(PetscObjectChangeTypeName((PetscObject)field, type));
    127:   field->ops->create = r;
    128:   PetscCall((*r)(field));
    129:   PetscFunctionReturn(PETSC_SUCCESS);
    130: }

    132: /*@
    133:   DMFieldGetType - Gets the `DMFieldType` name (as a string) from the `DMField`.

    135:   Not Collective

    137:   Input Parameter:
    138: . field - The `DMField` context

    140:   Output Parameter:
    141: . type - The `DMFieldType` name

    143:   Level: advanced

    145: .seealso: `DMField`, `DMFieldSetType()`, `DMFieldType`
    146: @*/
    147: PetscErrorCode DMFieldGetType(DMField field, DMFieldType *type)
    148: {
    149:   PetscFunctionBegin;
    151:   PetscAssertPointer(type, 2);
    152:   PetscCall(DMFieldRegisterAll());
    153:   *type = ((PetscObject)field)->type_name;
    154:   PetscFunctionReturn(PETSC_SUCCESS);
    155: }

    157: /*@
    158:   DMFieldGetNumComponents - Returns the number of components in the field

    160:   Not Collective

    162:   Input Parameter:
    163: . field - The `DMField` object

    165:   Output Parameter:
    166: . nc - The number of field components

    168:   Level: intermediate

    170: .seealso: `DMField`, `DMFieldEvaluate()`
    171: @*/
    172: PetscErrorCode DMFieldGetNumComponents(DMField field, PetscInt *nc)
    173: {
    174:   PetscFunctionBegin;
    176:   PetscAssertPointer(nc, 2);
    177:   *nc = field->numComponents;
    178:   PetscFunctionReturn(PETSC_SUCCESS);
    179: }

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

    184:   Not Collective

    186:   Input Parameter:
    187: . field - The `DMField` object

    189:   Output Parameter:
    190: . dm - The `DM` object

    192:   Level: intermediate

    194: .seealso: `DMField`, `DM`, `DMFieldEvaluate()`
    195: @*/
    196: PetscErrorCode DMFieldGetDM(DMField field, DM *dm)
    197: {
    198:   PetscFunctionBegin;
    200:   PetscAssertPointer(dm, 2);
    201:   *dm = field->dm;
    202:   PetscFunctionReturn(PETSC_SUCCESS);
    203: }

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

    208:   Collective

    210:   Input Parameters:
    211: + field    - The `DMField` object
    212: . points   - The points at which to evaluate the field.  Should have size d x n,
    213:            where d is the coordinate dimension of the manifold and n is the number
    214:            of points
    215: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
    216:              If the field is complex and datatype is `PETSC_REAL`, the real part of the
    217:              field is returned.

    219:   Output Parameters:
    220: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
    221:       If B is not NULL, the values of the field are written in this array, varying first by component,
    222:       then by point.
    223: . D - pointer to data of size d * c * n * sizeof(datatype).
    224:       If D is not NULL, the values of the field's spatial derivatives are written in this array,
    225:       varying first by the partial derivative component, then by field component, then by point.
    226: - H - pointer to data of size d * d * c * n * sizeof(datatype).
    227:       If H is not NULL, the values of the field's second spatial derivatives are written in this array,
    228:       varying first by the second partial derivative component, then by field component, then by point.

    230:   Level: intermediate

    232: .seealso: `DMField`, `DMFieldGetDM()`, `DMFieldGetNumComponents()`, `DMFieldEvaluateFE()`, `DMFieldEvaluateFV()`, `PetscDataType`
    233: @*/
    234: PetscErrorCode DMFieldEvaluate(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
    235: {
    236:   PetscFunctionBegin;
    239:   if (B) PetscAssertPointer(B, 4);
    240:   if (D) PetscAssertPointer(D, 5);
    241:   if (H) PetscAssertPointer(H, 6);
    242:   PetscUseTypeMethod(field, evaluate, points, datatype, B, D, H);
    243:   PetscFunctionReturn(PETSC_SUCCESS);
    244: }

    246: /*@
    247:   DMFieldEvaluateFE - Evaluate the field and its derivatives on a set of points mapped from
    248:   quadrature points on a reference point.  The derivatives are taken with respect to the
    249:   reference coordinates.

    251:   Not Collective

    253:   Input Parameters:
    254: + field    - The `DMField` object
    255: . cellIS   - Index set for cells on which to evaluate the field
    256: . points   - The quadature containing the points in the reference cell at which to evaluate the field.
    257: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
    258:              If the field is complex and datatype is `PETSC_REAL`, the real part of the
    259:              field is returned.

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

    272:   Level: intermediate

    274: .seealso: `DMField`, `DM`, `DMFieldGetNumComponents()`, `DMFieldEvaluate()`, `DMFieldEvaluateFV()`
    275: @*/
    276: PetscErrorCode DMFieldEvaluateFE(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
    277: {
    278:   PetscFunctionBegin;
    282:   if (B) PetscAssertPointer(B, 5);
    283:   if (D) PetscAssertPointer(D, 6);
    284:   if (H) PetscAssertPointer(H, 7);
    285:   PetscUseTypeMethod(field, evaluateFE, cellIS, points, datatype, B, D, H);
    286:   PetscFunctionReturn(PETSC_SUCCESS);
    287: }

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

    292:   Not Collective

    294:   Input Parameters:
    295: + field    - The `DMField` object
    296: . cellIS   - Index set for cells on which to evaluate the field
    297: - datatype - The PetscDataType of the output arrays: either `PETSC_REAL` or `PETSC_SCALAR`.
    298:              If the field is complex and datatype is `PETSC_REAL`, the real part of the
    299:              field is returned.

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

    312:   Level: intermediate

    314: .seealso: `DMField`, `IS`, `DMFieldGetNumComponents()`, `DMFieldEvaluate()`, `DMFieldEvaluateFE()`, `PetscDataType`
    315: @*/
    316: PetscErrorCode DMFieldEvaluateFV(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
    317: {
    318:   PetscFunctionBegin;
    321:   if (B) PetscAssertPointer(B, 4);
    322:   if (D) PetscAssertPointer(D, 5);
    323:   if (H) PetscAssertPointer(H, 6);
    324:   PetscUseTypeMethod(field, evaluateFV, cellIS, datatype, B, D, H);
    325:   PetscFunctionReturn(PETSC_SUCCESS);
    326: }

    328: /*@
    329:   DMFieldGetDegree - Get the polynomial degree of a field when pulled back onto the
    330:   reference element

    332:   Not Collective

    334:   Input Parameters:
    335: + field  - the `DMField` object
    336: - cellIS - the index set of points over which we want know the invariance

    338:   Output Parameters:
    339: + minDegree - the degree of the largest polynomial space contained in the field on each element
    340: - maxDegree - the largest degree of the smallest polynomial space containing the field on any element

    342:   Level: intermediate

    344: .seealso: `DMField`, `IS`, `DMFieldEvaluateFE()`
    345: @*/
    346: PetscErrorCode DMFieldGetDegree(DMField field, IS cellIS, PetscInt *minDegree, PetscInt *maxDegree)
    347: {
    348:   PetscFunctionBegin;
    351:   if (minDegree) PetscAssertPointer(minDegree, 3);
    352:   if (maxDegree) PetscAssertPointer(maxDegree, 4);

    354:   if (minDegree) *minDegree = -1;
    355:   if (maxDegree) *maxDegree = PETSC_INT_MAX;

    357:   PetscTryTypeMethod(field, getDegree, cellIS, minDegree, maxDegree);
    358:   PetscFunctionReturn(PETSC_SUCCESS);
    359: }

    361: /*@
    362:   DMFieldCreateDefaultQuadrature - Creates a quadrature sufficient to integrate the field on the selected
    363:   points via pullback onto the reference element

    365:   Not Collective

    367:   Input Parameters:
    368: + field   - the `DMField` object
    369: - pointIS - the index set of points over which we wish to integrate the field

    371:   Output Parameter:
    372: . quad - a `PetscQuadrature` object

    374:   Level: developer

    376: .seealso: `DMField`, `PetscQuadrature`, `IS`, `DMFieldEvaluteFE()`, `DMFieldGetDegree()`
    377: @*/
    378: PetscErrorCode DMFieldCreateDefaultQuadrature(DMField field, IS pointIS, PetscQuadrature *quad)
    379: {
    380:   PetscFunctionBegin;
    383:   PetscAssertPointer(quad, 3);

    385:   *quad = NULL;
    386:   PetscTryTypeMethod(field, createDefaultQuadrature, pointIS, quad);
    387:   PetscFunctionReturn(PETSC_SUCCESS);
    388: }

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

    393:   Not Collective

    395:   Input Parameters:
    396: + field    - the `DMField` object
    397: . pointIS  - the index set of points over which we wish to integrate the field
    398: . quad     - the quadrature points at which to evaluate the geometric factors
    399: - faceData - whether additional data for facets (the normal vectors and adjacent cells) should
    400:   be calculated

    402:   Output Parameter:
    403: . geom - the geometric factors

    405:   Level: developer

    407: .seealso: `DMField`, `PetscQuadrature`, `IS`, `PetscFEGeom`, `DMFieldEvaluateFE()`, `DMFieldCreateDefaulteQuadrature()`, `DMFieldGetDegree()`
    408: @*/
    409: PetscErrorCode DMFieldCreateFEGeom(DMField field, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
    410: {
    411:   PetscInt     dim, dE;
    412:   PetscInt     nPoints;
    413:   PetscInt     maxDegree;
    414:   PetscFEGeom *g;

    416:   PetscFunctionBegin;
    420:   PetscCall(ISGetLocalSize(pointIS, &nPoints));
    421:   dE = field->numComponents;
    422:   PetscCall(PetscFEGeomCreate(quad, nPoints, dE, faceData, &g));
    423:   PetscCall(DMFieldEvaluateFE(field, pointIS, quad, PETSC_REAL, g->v, g->J, NULL));
    424:   dim = g->dim;
    425:   if (dE > dim) {
    426:     /* space out J and make square Jacobians */
    427:     PetscInt i, j, k, N = g->numPoints * g->numCells;

    429:     for (i = N - 1; i >= 0; i--) {
    430:       PetscReal J[9] = {0};

    432:       for (j = 0; j < dE; j++) {
    433:         for (k = 0; k < dim; k++) J[j * dE + k] = g->J[i * dE * dim + j * dim + k];
    434:       }
    435:       switch (dim) {
    436:       case 0:
    437:         for (j = 0; j < dE; j++) {
    438:           for (k = 0; k < dE; k++) J[j * dE + k] = (j == k) ? 1. : 0.;
    439:         }
    440:         break;
    441:       case 1:
    442:         if (dE == 2) {
    443:           PetscReal norm = PetscSqrtReal(J[0] * J[0] + J[2] * J[2]);

    445:           J[1] = -J[2] / norm;
    446:           J[3] = J[0] / norm;
    447:         } else {
    448:           PetscReal inorm = 1. / PetscSqrtReal(J[0] * J[0] + J[3] * J[3] + J[6] * J[6]);
    449:           PetscReal x     = J[0] * inorm;
    450:           PetscReal y     = J[3] * inorm;
    451:           PetscReal z     = J[6] * inorm;

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

    456:             J[1] = -y;
    457:             J[2] = -z;
    458:             J[4] = 1. - y * y * inv1pX;
    459:             J[5] = -y * z * inv1pX;
    460:             J[7] = -y * z * inv1pX;
    461:             J[8] = 1. - z * z * inv1pX;
    462:           } else {
    463:             PetscReal inv1mX = 1. / (1. - x);

    465:             J[1] = z;
    466:             J[2] = y;
    467:             J[4] = -y * z * inv1mX;
    468:             J[5] = 1. - y * y * inv1mX;
    469:             J[7] = 1. - z * z * inv1mX;
    470:             J[8] = -y * z * inv1mX;
    471:           }
    472:         }
    473:         break;
    474:       case 2: {
    475:         PetscReal inorm;

    477:         J[2] = J[3] * J[7] - J[6] * J[4];
    478:         J[5] = J[6] * J[1] - J[0] * J[7];
    479:         J[8] = J[0] * J[4] - J[3] * J[1];

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

    483:         J[2] *= inorm;
    484:         J[5] *= inorm;
    485:         J[8] *= inorm;
    486:       } break;
    487:       }
    488:       for (j = 0; j < dE * dE; j++) g->J[i * dE * dE + j] = J[j];
    489:     }
    490:   }
    491:   PetscCall(PetscFEGeomComplete(g));
    492:   PetscCall(DMFieldGetDegree(field, pointIS, NULL, &maxDegree));
    493:   g->isAffine = (maxDegree <= 1) ? PETSC_TRUE : PETSC_FALSE;
    494:   if (faceData) PetscUseTypeMethod(field, computeFaceData, pointIS, quad, g);
    495:   *geom = g;
    496:   PetscFunctionReturn(PETSC_SUCCESS);
    497: }