Actual source code: dmfieldshell.c

petsc-3.12.2 2019-11-22
Report Typos and Errors
  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: }