Actual source code: dmfieldshell.c

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

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

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

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

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

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

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

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

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

 58:           PetscInt i, j;

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

 71:       for (p = 0; p < numPoints * Nq; p++) {
 72:         for (q = 0; q < Nc; q++) {
 73:           PetscReal d[3];

 75:           PetscInt i, j;

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

 91:       for (p = 0; p < numPoints * Nq; p++) {
 92:         for (q = 0; q < Nc; q++) {
 93:           PetscScalar d[3][3];

 95:           PetscInt i, j, k, l;

 97:           for (i = 0; i < dimC; i++)
 98:             for (j = 0; j < dimC; j++) d[i][j] = sH[((p * Nc + q) * dimC + i) * dimC + j];
 99:           for (i = 0; i < dimC; i++)
100:             for (j = 0; j < dimC; j++) sH[((p * Nc + q) * dimC + i) * dimC + j] = 0.;
101:           for (i = 0; i < dimC; i++) {
102:             for (j = 0; j < dimC; j++) {
103:               for (k = 0; k < dimC; k++) {
104:                 for (l = 0; l < dimC; l++) 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];
105:               }
106:             }
107:           }
108:         }
109:       }
110:     } else {
111:       PetscReal *rH = (PetscReal *)H;
112:       PetscInt   q;

114:       for (p = 0; p < numPoints * Nq; p++) {
115:         for (q = 0; q < Nc; q++) {
116:           PetscReal d[3][3];

118:           PetscInt i, j, k, l;

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

141: PetscErrorCode DMFieldShellEvaluateFVDefault(DMField field, IS pointIS, PetscDataType type, void *B, void *D, void *H)
142: {
143:   DM              dm = field->dm;
144:   DMField         coordField;
145:   PetscFEGeom    *geom;
146:   Vec             pushforward;
147:   PetscInt        dimC, dim, numPoints, Nq, p;
148:   PetscScalar    *pfArray;
149:   PetscQuadrature quad;
150:   MPI_Comm        comm;

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

173: PetscErrorCode DMFieldShellSetDestroy(DMField field, PetscErrorCode (*destroy)(DMField))
174: {
175:   DMField_Shell *shell = (DMField_Shell *)field->data;

178:   shell->destroy = destroy;
179:   return 0;
180: }

182: PetscErrorCode DMFieldShellSetEvaluate(DMField field, PetscErrorCode (*evaluate)(DMField, Vec, PetscDataType, void *, void *, void *))
183: {
185:   field->ops->evaluate = evaluate;
186:   return 0;
187: }

189: PetscErrorCode DMFieldShellSetEvaluateFE(DMField field, PetscErrorCode (*evaluateFE)(DMField, IS, PetscQuadrature, PetscDataType, void *, void *, void *))
190: {
192:   field->ops->evaluateFE = evaluateFE;
193:   return 0;
194: }

196: PetscErrorCode DMFieldShellSetEvaluateFV(DMField field, PetscErrorCode (*evaluateFV)(DMField, IS, PetscDataType, void *, void *, void *))
197: {
199:   field->ops->evaluateFV = evaluateFV;
200:   return 0;
201: }

203: PetscErrorCode DMFieldShellSetGetDegree(DMField field, PetscErrorCode (*getDegree)(DMField, IS, PetscInt *, PetscInt *))
204: {
206:   field->ops->getDegree = getDegree;
207:   return 0;
208: }

210: PetscErrorCode DMFieldShellSetCreateDefaultQuadrature(DMField field, PetscErrorCode (*createDefaultQuadrature)(DMField, IS, PetscQuadrature *))
211: {
213:   field->ops->createDefaultQuadrature = createDefaultQuadrature;
214:   return 0;
215: }

217: static PetscErrorCode DMFieldInitialize_Shell(DMField field)
218: {
219:   field->ops->destroy                 = DMFieldDestroy_Shell;
220:   field->ops->evaluate                = NULL;
221:   field->ops->evaluateFE              = DMFieldShellEvaluateFEDefault;
222:   field->ops->evaluateFV              = DMFieldShellEvaluateFVDefault;
223:   field->ops->getDegree               = NULL;
224:   field->ops->createDefaultQuadrature = NULL;
225:   field->ops->view                    = NULL;
226:   return 0;
227: }

229: PETSC_INTERN PetscErrorCode DMFieldCreate_Shell(DMField field)
230: {
231:   DMField_Shell *shell;

233:   PetscNew(&shell);
234:   field->data = shell;
235:   DMFieldInitialize_Shell(field);
236:   return 0;
237: }

239: PetscErrorCode DMFieldCreateShell(DM dm, PetscInt numComponents, DMFieldContinuity continuity, void *ctx, DMField *field)
240: {
241:   DMField        b;
242:   DMField_Shell *shell;

247:   DMFieldCreate(dm, numComponents, continuity, &b);
248:   DMFieldSetType(b, DMFIELDSHELL);
249:   shell      = (DMField_Shell *)b->data;
250:   shell->ctx = ctx;
251:   *field     = b;
252:   return 0;
253: }