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[] = {
  6:   "VERTEX",
  7:   "EDGE",
  8:   "FACET",
  9:   "CELL",
 10:   NULL
 11: };

 13: PETSC_INTERN PetscErrorCode DMFieldCreate(DM dm,PetscInt numComponents,DMFieldContinuity continuity,DMField *field)
 14: {
 15:   DMField        b;

 19:   DMFieldInitializePackage();

 21:   PetscHeaderCreate(b,DMFIELD_CLASSID,"DMField","Field over DM","DM",PetscObjectComm((PetscObject)dm),DMFieldDestroy,DMFieldView);
 22:   PetscObjectReference((PetscObject)dm);
 23:   b->dm = dm;
 24:   b->continuity = continuity;
 25:   b->numComponents = numComponents;
 26:   *field = b;
 27:   return 0;
 28: }

 30: /*@
 31:    DMFieldDestroy - destroy a DMField

 33:    Collective

 35:    Input Parameter:
 36: .  field - address of DMField

 38:    Level: advanced

 40: .seealso: DMFieldCreate()
 41: @*/
 42: PetscErrorCode DMFieldDestroy(DMField *field)
 43: {
 44:   if (!*field) return 0;
 46:   if (--((PetscObject)(*field))->refct > 0) {*field = NULL; return 0;}
 47:   if ((*field)->ops->destroy) (*(*field)->ops->destroy)(*field);
 48:   DMDestroy(&((*field)->dm));
 49:   PetscHeaderDestroy(field);
 50:   return 0;
 51: }

 53: /*@C
 54:    DMFieldView - view a DMField

 56:    Collective

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

 62:    Level: advanced

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

 71:   if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)field),&viewer);
 74:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 75:   if (iascii) {
 76:     PetscObjectPrintClassNamePrefixType((PetscObject)field,viewer);
 77:     PetscViewerASCIIPushTab(viewer);
 78:     PetscViewerASCIIPrintf(viewer,"%D components\n",field->numComponents);
 79:     PetscViewerASCIIPrintf(viewer,"%s continuity\n",DMFieldContinuities[field->continuity]);
 80:     PetscViewerPushFormat(viewer,PETSC_VIEWER_DEFAULT);
 81:     DMView(field->dm,viewer);
 82:     PetscViewerPopFormat(viewer);
 83:   }
 84:   if (field->ops->view) (*field->ops->view)(field,viewer);
 85:   if (iascii) {
 86:     PetscViewerASCIIPopTab(viewer);
 87:   }
 88:   return 0;
 89: }

 91: /*@C
 92:    DMFieldSetType - set the DMField implementation

 94:    Collective on field

 96:    Input Parameters:
 97: +  field - the DMField context
 98: -  type - a known method

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:   Level: advanced

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


118:   PetscObjectTypeCompare((PetscObject)field,type,&match);
119:   if (match) return 0;

121:   PetscFunctionListFind(DMFieldList,type,&r);
123:   /* Destroy the previous private DMField context */
124:   if (field->ops->destroy) (*(field)->ops->destroy)(field);

126:   PetscMemzero(field->ops,sizeof(*field->ops));
127:   PetscObjectChangeTypeName((PetscObject)field,type);
128:   field->ops->create = r;
129:   (*r)(field);
130:   return 0;
131: }

133: /*@C
134:   DMFieldGetType - Gets the DMField type name (as a string) from the DMField.

136:   Not Collective

138:   Input Parameter:
139: . field  - The DMField context

141:   Output Parameter:
142: . type - The DMField type name

144:   Level: advanced

146: .seealso: DMFieldSetType()
147: @*/
148: PetscErrorCode  DMFieldGetType(DMField field, DMFieldType *type)
149: {
152:   DMFieldRegisterAll();
153:   *type = ((PetscObject)field)->type_name;
154:   return 0;
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: DMFieldEvaluate()
171: @*/
172: PetscErrorCode DMFieldGetNumComponents(DMField field, PetscInt *nc)
173: {
176:   *nc = field->numComponents;
177:   return 0;
178: }

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

183:   Not collective

185:   Input Parameter:
186: . field - The DMField object

188:   Output Parameter:
189: . dm - The DM object

191:   Level: intermediate

193: .seealso: DMFieldEvaluate()
194: @*/
195: PetscErrorCode DMFieldGetDM(DMField field, DM *dm)
196: {
199:   *dm = field->dm;
200:   return 0;
201: }

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

206:   Collective on points

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

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

228:   Level: intermediate

230: .seealso: DMFieldGetDM(), DMFieldGetNumComponents(), DMFieldEvaluateFE(), DMFieldEvaluateFV()
231: @*/
232: PetscErrorCode DMFieldEvaluate(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
233: {
239:   if (field->ops->evaluate) {
240:     (*field->ops->evaluate) (field, points, datatype, B, D, H);
241:   } else SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
242:   return 0;
243: }

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

250:   Not collective

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

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

271:   Level: intermediate

273: .seealso: DMFieldGetNumComponents(), DMFieldEvaluate(), DMFieldEvaluateFV()
274: @*/
275: PetscErrorCode DMFieldEvaluateFE(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
276: {
283:   if (field->ops->evaluateFE) {
284:     (*field->ops->evaluateFE) (field, cellIS, points, datatype, B, D, H);
285:   } else SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
286:   return 0;
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: DMFieldGetNumComponents(), DMFieldEvaluate(), DMFieldEvaluateFE()
315: @*/
316: PetscErrorCode DMFieldEvaluateFV(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
317: {
323:   if (field->ops->evaluateFV) {
324:     (*field->ops->evaluateFV) (field, cellIS, datatype, B, D, H);
325:   } else SETERRQ(PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
326:   return 0;
327: }

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

333:   Not collective

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

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

343:   Level: intermediate

345: .seealso: DMFieldEvaluateFE()
346: @*/
347: PetscErrorCode DMFieldGetDegree(DMField field, IS cellIS, PetscInt *minDegree, PetscInt *maxDegree)
348: {

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

357:   if (field->ops->getDegree) {
358:     (*field->ops->getDegree) (field,cellIS,minDegree,maxDegree);
359:   }
360:   return 0;
361: }

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

367:   Not collective

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

373:   Output Parameter:
374: . quad - a PetscQuadrature object

376:   Level: developer

378: .seealso: DMFieldEvaluteFE(), DMFieldGetDegree()
379: @*/
380: PetscErrorCode DMFieldCreateDefaultQuadrature(DMField field, IS pointIS, PetscQuadrature *quad)
381: {

386:   *quad = NULL;
387:   if (field->ops->createDefaultQuadrature) {
388:     (*field->ops->createDefaultQuadrature)(field, pointIS, quad);
389:   }
390:   return 0;
391: }

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

396:   Not collective

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

405:   Output Parameter:
406: . geom - the geometric factors

408:   Level: developer

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

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

431:     for (i = N-1; i >= 0; i--) {
432:       PetscReal   J[9];

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

451:           J[1] = -J[2] / norm;
452:           J[3] =  J[0] / norm;
453:         } else {
454:           PetscReal inorm = 1./PetscSqrtReal(J[0] * J[0] + J[3] * J[3] + J[6] * J[6]);
455:           PetscReal x = J[0] * inorm;
456:           PetscReal y = J[3] * inorm;
457:           PetscReal z = J[6] * inorm;

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

462:             J[1] = -y;              J[2] = -z;
463:             J[4] = 1. - y*y*inv1pX; J[5] =     -y*z*inv1pX;
464:             J[7] =     -y*z*inv1pX; J[8] = 1. - z*z*inv1pX;
465:           } else {
466:             PetscReal inv1mX   = 1./ (1. - x);

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

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

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

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