Actual source code: dmfieldshell.c

      1: #include <petsc/private/dmfieldimpl.h>

      3: typedef struct _n_DMField_Shell
      4: {
      5:   void *ctx;
      6:   PetscErrorCode (*destroy) (DMField);
      7: }
      8: DMField_Shell;

     10: PetscErrorCode DMFieldShellGetContext(DMField field, void *ctx)
     11: {
     12:   PetscBool      flg;

     18:   PetscObjectTypeCompare((PetscObject)field,DMFIELDSHELL,&flg);
     19:   if (flg) *(void**)ctx = ((DMField_Shell*)(field->data))->ctx;
     20:   else SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Cannot get context from non-shell shield");
     21:   return(0);
     22: }

     24: static PetscErrorCode DMFieldDestroy_Shell(DMField field)
     25: {
     26:   DMField_Shell *shell = (DMField_Shell *) field->data;

     30:   if (shell->destroy) {(*(shell->destroy)) (field);}
     31:   PetscFree(field->data);
     32:   return(0);
     33: }

     35: PetscErrorCode DMFieldShellEvaluateFEDefault(DMField field, IS pointIS, PetscQuadrature quad, PetscDataType type, void *B, void *D, void *H)
     36: {
     37:   DM              dm = field->dm;
     38:   DMField         coordField;
     39:   PetscFEGeom    *geom;
     40:   Vec             pushforward;
     41:   PetscInt        dimC, dim, numPoints, Nq, p, Nc;
     42:   PetscScalar    *pfArray;
     43:   PetscErrorCode  ierr;

     46:   Nc   = field->numComponents;
     47:   DMGetCoordinateField(dm, &coordField);
     48:   DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom);
     49:   DMGetCoordinateDim(dm, &dimC);
     50:   PetscQuadratureGetData(quad, &dim, NULL, &Nq, NULL, NULL);
     51:   ISGetLocalSize(pointIS, &numPoints);
     52:   PetscMalloc1(dimC * Nq * numPoints, &pfArray);
     53:   for (p = 0; p < numPoints * dimC * Nq; p++) pfArray[p] = (PetscScalar) geom->v[p];
     54:   VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * Nq * numPoints, PETSC_DETERMINE, pfArray, &pushforward);
     55:   DMFieldEvaluate(field, pushforward, type, B, D, H);
     56:   /* TODO: handle covariant/contravariant pullbacks */
     57:   if (D) {
     58:     if (type == PETSC_SCALAR) {
     59:       PetscScalar *sD = (PetscScalar *) D;
     60:       PetscInt q;

     62:       for (p = 0; p < numPoints * Nq; p++) {
     63:         for (q = 0; q < Nc; q++) {
     64:           PetscScalar d[3];

     66:           PetscInt i, j;

     68:           for (i = 0; i < dimC; i++) d[i] = sD[(p * Nc + q) * dimC + i];
     69:           for (i = 0; i < dimC; i++) sD[(p * Nc + q) * dimC + i] = 0.;
     70:           for (i = 0; i < dimC; i++) {
     71:             for (j = 0; j < dimC; j++) {
     72:               sD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
     73:             }
     74:           }
     75:         }
     76:       }
     77:     } else {
     78:       PetscReal *rD = (PetscReal *) D;
     79:       PetscInt q;

     81:       for (p = 0; p < numPoints * Nq; p++) {
     82:         for (q = 0; q < Nc; q++) {
     83:           PetscReal d[3];

     85:           PetscInt i, j;

     87:           for (i = 0; i < dimC; i++) d[i] = rD[(p * Nc + q) * dimC + i];
     88:           for (i = 0; i < dimC; i++) rD[(p * Nc + q) * dimC + i] = 0.;
     89:           for (i = 0; i < dimC; i++) {
     90:             for (j = 0; j < dimC; j++) {
     91:               rD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
     92:             }
     93:           }
     94:         }
     95:       }
     96:     }
     97:   }
     98:   if (H) {
     99:     if (type == PETSC_SCALAR) {
    100:       PetscScalar *sH = (PetscScalar *) H;
    101:       PetscInt q;

    103:       for (p = 0; p < numPoints * Nq; p++) {
    104:         for (q = 0; q < Nc; q++) {
    105:           PetscScalar d[3][3];

    107:           PetscInt i, j, k, l;

    109:           for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) d[i][j] = sH[((p * Nc + q) * dimC + i) * dimC + j];
    110:           for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) sH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
    111:           for (i = 0; i < dimC; i++) {
    112:             for (j = 0; j < dimC; j++) {
    113:               for (k = 0; k < dimC; k++) {
    114:                 for (l = 0; l < dimC; l++) {
    115:                   sH[((p * Nc + q) * dimC + i) * dimC + j] += geom->J[(p * dimC + k) * dimC + i] * geom->J[(p * dimC + l) * dimC + j] * d[k][l];
    116:                 }
    117:               }
    118:             }
    119:           }
    120:         }
    121:       }
    122:     } else {
    123:       PetscReal *rH = (PetscReal *) H;
    124:       PetscInt q;

    126:       for (p = 0; p < numPoints * Nq; p++) {
    127:         for (q = 0; q < Nc; q++) {
    128:           PetscReal d[3][3];

    130:           PetscInt i, j, k, l;

    132:           for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) d[i][j] = rH[((p * Nc + q) * dimC + i) * dimC + j];
    133:           for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) rH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
    134:           for (i = 0; i < dimC; i++) {
    135:             for (j = 0; j < dimC; j++) {
    136:               for (k = 0; k < dimC; k++) {
    137:                 for (l = 0; l < dimC; l++) {
    138:                   rH[((p * Nc + q) * dimC + i) * dimC + j] += geom->J[(p * dimC + k) * dimC + i] * geom->J[(p * dimC + l) * dimC + j] * d[k][l];
    139:                 }
    140:               }
    141:             }
    142:           }
    143:         }
    144:       }
    145:     }
    146:   }
    147:   VecDestroy(&pushforward);
    148:   PetscFree(pfArray);
    149:   PetscFEGeomDestroy(&geom);
    150:   return(0);
    151: }

    153: PetscErrorCode DMFieldShellEvaluateFVDefault(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
    154: {
    155:   DM              dm = field->dm;
    156:   DMField         coordField;
    157:   PetscFEGeom    *geom;
    158:   Vec             pushforward;
    159:   PetscInt        dimC, dim, numPoints, Nq, p;
    160:   PetscScalar    *pfArray;
    161:   PetscQuadrature quad;
    162:   PetscErrorCode  ierr;

    165:   DMGetCoordinateField(dm, &coordField);
    166:   DMFieldCreateDefaultQuadrature(coordField, pointIS, &quad);
    167:   if (!quad) SETERRQ(PetscObjectComm((PetscObject) pointIS), PETSC_ERR_ARG_WRONGSTATE, "coordinate field must have default quadrature for FV computation");
    168:   DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom);
    169:   DMGetCoordinateDim(dm, &dimC);
    170:   PetscQuadratureGetData(quad, &dim, NULL, &Nq, NULL, NULL);
    171:   if (Nq != 1) SETERRQ(PetscObjectComm((PetscObject) quad), PETSC_ERR_ARG_WRONGSTATE, "quadrature must have only one point");
    172:   ISGetLocalSize(pointIS, &numPoints);
    173:   PetscMalloc1(dimC * numPoints, &pfArray);
    174:   for (p = 0; p < numPoints * dimC; p++) pfArray[p] = (PetscScalar) geom->v[p];
    175:   VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * numPoints, PETSC_DETERMINE, pfArray, &pushforward);
    176:   DMFieldEvaluate(field, pushforward, type, B, D, H);
    177:   PetscQuadratureDestroy(&quad);
    178:   VecDestroy(&pushforward);
    179:   PetscFree(pfArray);
    180:   PetscFEGeomDestroy(&geom);
    181:   return(0);
    182: }

    184: PetscErrorCode DMFieldShellSetDestroy(DMField field, PetscErrorCode (*destroy)(DMField))
    185: {
    186:   DMField_Shell *shell = (DMField_Shell *) field->data;

    190:   shell->destroy = destroy;
    191:   return(0);
    192: }

    194: PetscErrorCode DMFieldShellSetEvaluate(DMField field, PetscErrorCode (*evaluate)(DMField,Vec,PetscDataType,void*,void*,void*))
    195: {
    198:   field->ops->evaluate = evaluate;
    199:   return(0);
    200: }

    202: PetscErrorCode DMFieldShellSetEvaluateFE(DMField field, PetscErrorCode (*evaluateFE)(DMField,IS,PetscQuadrature,PetscDataType,void*,void*,void*))
    203: {
    206:   field->ops->evaluateFE = evaluateFE;
    207:   return(0);
    208: }

    210: PetscErrorCode DMFieldShellSetEvaluateFV(DMField field, PetscErrorCode (*evaluateFV)(DMField,IS,PetscDataType,void*,void*,void*))
    211: {
    214:   field->ops->evaluateFV = evaluateFV;
    215:   return(0);
    216: }

    218: PetscErrorCode DMFieldShellSetGetDegree(DMField field, PetscErrorCode (*getDegree)(DMField,IS,PetscInt*,PetscInt*))
    219: {
    222:   field->ops->getDegree = getDegree;
    223:   return(0);
    224: }

    226: PetscErrorCode DMFieldShellSetCreateDefaultQuadrature(DMField field, PetscErrorCode (*createDefaultQuadrature)(DMField,IS,PetscQuadrature*))
    227: {
    230:   field->ops->createDefaultQuadrature = createDefaultQuadrature;
    231:   return(0);
    232: }

    234: static PetscErrorCode DMFieldInitialize_Shell(DMField field)
    235: {
    237:   field->ops->destroy                 = DMFieldDestroy_Shell;
    238:   field->ops->evaluate                = NULL;
    239:   field->ops->evaluateFE              = DMFieldShellEvaluateFEDefault;
    240:   field->ops->evaluateFV              = DMFieldShellEvaluateFVDefault;
    241:   field->ops->getDegree               = NULL;
    242:   field->ops->createDefaultQuadrature = NULL;
    243:   field->ops->view                    = NULL;
    244:   return(0);
    245: }

    247: PETSC_INTERN PetscErrorCode DMFieldCreate_Shell(DMField field)
    248: {
    249:   DMField_Shell *shell;

    253:   PetscNewLog(field,&shell);
    254:   field->data = shell;
    255:   DMFieldInitialize_Shell(field);
    256:   return(0);
    257: }

    259: PetscErrorCode DMFieldCreateShell(DM dm, PetscInt numComponents, DMFieldContinuity continuity, void *ctx, DMField *field)
    260: {
    261:   DMField        b;
    262:   DMField_Shell  *shell;

    269:   DMFieldCreate(dm,numComponents,continuity,&b);
    270:   DMFieldSetType(b,DMFIELDSHELL);
    271:   shell = (DMField_Shell *) b->data;
    272:   shell->ctx = ctx;
    273:   *field = b;
    274:   return(0);
    275: }