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;

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

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

 26:   if (shell->destroy) (*(shell->destroy)) (field);
 27:   PetscFree(field->data);
 28:   return 0;
 29: }

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

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

 56:       for (p = 0; p < numPoints * Nq; p++) {
 57:         for (q = 0; q < Nc; q++) {
 58:           PetscScalar d[3];

 60:           PetscInt i, j;

 62:           for (i = 0; i < dimC; i++) d[i] = sD[(p * Nc + q) * dimC + i];
 63:           for (i = 0; i < dimC; i++) sD[(p * Nc + q) * dimC + i] = 0.;
 64:           for (i = 0; i < dimC; i++) {
 65:             for (j = 0; j < dimC; j++) {
 66:               sD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
 67:             }
 68:           }
 69:         }
 70:       }
 71:     } else {
 72:       PetscReal *rD = (PetscReal *) D;
 73:       PetscInt q;

 75:       for (p = 0; p < numPoints * Nq; p++) {
 76:         for (q = 0; q < Nc; q++) {
 77:           PetscReal d[3];

 79:           PetscInt i, j;

 81:           for (i = 0; i < dimC; i++) d[i] = rD[(p * Nc + q) * dimC + i];
 82:           for (i = 0; i < dimC; i++) rD[(p * Nc + q) * dimC + i] = 0.;
 83:           for (i = 0; i < dimC; i++) {
 84:             for (j = 0; j < dimC; j++) {
 85:               rD[(p * Nc + q) * dimC + i] += geom->J[(p * dimC + j) * dimC + i] * d[j];
 86:             }
 87:           }
 88:         }
 89:       }
 90:     }
 91:   }
 92:   if (H) {
 93:     if (type == PETSC_SCALAR) {
 94:       PetscScalar *sH = (PetscScalar *) H;
 95:       PetscInt q;

 97:       for (p = 0; p < numPoints * Nq; p++) {
 98:         for (q = 0; q < Nc; q++) {
 99:           PetscScalar d[3][3];

101:           PetscInt i, j, k, l;

103:           for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) d[i][j] = sH[((p * Nc + q) * dimC + i) * dimC + j];
104:           for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) sH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
105:           for (i = 0; i < dimC; i++) {
106:             for (j = 0; j < dimC; j++) {
107:               for (k = 0; k < dimC; k++) {
108:                 for (l = 0; l < dimC; l++) {
109:                   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];
110:                 }
111:               }
112:             }
113:           }
114:         }
115:       }
116:     } else {
117:       PetscReal *rH = (PetscReal *) H;
118:       PetscInt q;

120:       for (p = 0; p < numPoints * Nq; p++) {
121:         for (q = 0; q < Nc; q++) {
122:           PetscReal d[3][3];

124:           PetscInt i, j, k, l;

126:           for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) d[i][j] = rH[((p * Nc + q) * dimC + i) * dimC + j];
127:           for (i = 0; i < dimC; i++) for (j = 0; j < dimC; j++) rH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
128:           for (i = 0; i < dimC; i++) {
129:             for (j = 0; j < dimC; j++) {
130:               for (k = 0; k < dimC; k++) {
131:                 for (l = 0; l < dimC; l++) {
132:                   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];
133:                 }
134:               }
135:             }
136:           }
137:         }
138:       }
139:     }
140:   }
141:   VecDestroy(&pushforward);
142:   PetscFree(pfArray);
143:   PetscFEGeomDestroy(&geom);
144:   return 0;
145: }

147: PetscErrorCode DMFieldShellEvaluateFVDefault(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
148: {
149:   DM              dm = field->dm;
150:   DMField         coordField;
151:   PetscFEGeom    *geom;
152:   Vec             pushforward;
153:   PetscInt        dimC, dim, numPoints, Nq, p;
154:   PetscScalar    *pfArray;
155:   PetscQuadrature quad;
156:   MPI_Comm        comm;

158:   PetscObjectGetComm((PetscObject) field, &comm);
159:   DMGetDimension(dm, &dim);
160:   DMGetCoordinateDim(dm, &dimC);
161:   DMGetCoordinateField(dm, &coordField);
162:   DMFieldGetFVQuadrature_Internal(coordField, pointIS, &quad);
164:   PetscQuadratureGetData(quad, NULL, NULL, &Nq, NULL, NULL);
166:   DMFieldCreateFEGeom(coordField, pointIS, quad, PETSC_FALSE, &geom);
167:   ISGetLocalSize(pointIS, &numPoints);
168:   PetscMalloc1(dimC * numPoints, &pfArray);
169:   for (p = 0; p < numPoints * dimC; p++) pfArray[p] = (PetscScalar) geom->v[p];
170:   VecCreateMPIWithArray(PetscObjectComm((PetscObject)pointIS), dimC, dimC * numPoints, PETSC_DETERMINE, pfArray, &pushforward);
171:   DMFieldEvaluate(field, pushforward, type, B, D, H);
172:   PetscQuadratureDestroy(&quad);
173:   VecDestroy(&pushforward);
174:   PetscFree(pfArray);
175:   PetscFEGeomDestroy(&geom);
176:   return 0;
177: }

179: PetscErrorCode DMFieldShellSetDestroy(DMField field, PetscErrorCode (*destroy)(DMField))
180: {
181:   DMField_Shell *shell = (DMField_Shell *) field->data;

184:   shell->destroy = destroy;
185:   return 0;
186: }

188: PetscErrorCode DMFieldShellSetEvaluate(DMField field, PetscErrorCode (*evaluate)(DMField,Vec,PetscDataType,void*,void*,void*))
189: {
191:   field->ops->evaluate = evaluate;
192:   return 0;
193: }

195: PetscErrorCode DMFieldShellSetEvaluateFE(DMField field, PetscErrorCode (*evaluateFE)(DMField,IS,PetscQuadrature,PetscDataType,void*,void*,void*))
196: {
198:   field->ops->evaluateFE = evaluateFE;
199:   return 0;
200: }

202: PetscErrorCode DMFieldShellSetEvaluateFV(DMField field, PetscErrorCode (*evaluateFV)(DMField,IS,PetscDataType,void*,void*,void*))
203: {
205:   field->ops->evaluateFV = evaluateFV;
206:   return 0;
207: }

209: PetscErrorCode DMFieldShellSetGetDegree(DMField field, PetscErrorCode (*getDegree)(DMField,IS,PetscInt*,PetscInt*))
210: {
212:   field->ops->getDegree = getDegree;
213:   return 0;
214: }

216: PetscErrorCode DMFieldShellSetCreateDefaultQuadrature(DMField field, PetscErrorCode (*createDefaultQuadrature)(DMField,IS,PetscQuadrature*))
217: {
219:   field->ops->createDefaultQuadrature = createDefaultQuadrature;
220:   return 0;
221: }

223: static PetscErrorCode DMFieldInitialize_Shell(DMField field)
224: {
225:   field->ops->destroy                 = DMFieldDestroy_Shell;
226:   field->ops->evaluate                = NULL;
227:   field->ops->evaluateFE              = DMFieldShellEvaluateFEDefault;
228:   field->ops->evaluateFV              = DMFieldShellEvaluateFVDefault;
229:   field->ops->getDegree               = NULL;
230:   field->ops->createDefaultQuadrature = NULL;
231:   field->ops->view                    = NULL;
232:   return 0;
233: }

235: PETSC_INTERN PetscErrorCode DMFieldCreate_Shell(DMField field)
236: {
237:   DMField_Shell *shell;

239:   PetscNewLog(field,&shell);
240:   field->data = shell;
241:   DMFieldInitialize_Shell(field);
242:   return 0;
243: }

245: PetscErrorCode DMFieldCreateShell(DM dm, PetscInt numComponents, DMFieldContinuity continuity, void *ctx, DMField *field)
246: {
247:   DMField        b;
248:   DMField_Shell  *shell;

253:   DMFieldCreate(dm,numComponents,continuity,&b);
254:   DMFieldSetType(b,DMFIELDSHELL);
255:   shell = (DMField_Shell *) b->data;
256:   shell->ctx = ctx;
257:   *field = b;
258:   return 0;
259: }