Actual source code: dmfield.c

petsc-master 2019-09-15
Report Typos and Errors
  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:   0
 11: };

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

 21:   DMFieldInitializePackage();

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

 32: /*@
 33:    DMFieldDestroy - destroy a DMField

 35:    Collective

 37:    Input Arguments:
 38: .  field - address of DMField

 40:    Level: advanced

 42: .seealso: DMFieldCreate()
 43: @*/
 44: PetscErrorCode DMFieldDestroy(DMField *field)
 45: {

 49:   if (!*field) return(0);
 51:   if (--((PetscObject)(*field))->refct > 0) {*field = 0; return(0);}
 52:   if ((*field)->ops->destroy) {(*(*field)->ops->destroy)(*field);}
 53:   DMDestroy(&((*field)->dm));
 54:   PetscHeaderDestroy(field);
 55:   return(0);
 56: }

 58: /*@C
 59:    DMFieldView - view a DMField

 61:    Collective

 63:    Input Arguments:
 64: +  field - DMField
 65: -  viewer - viewer to display field, for example PETSC_VIEWER_STDOUT_WORLD

 67:    Level: advanced

 69: .seealso: DMFieldCreate()
 70: @*/
 71: PetscErrorCode DMFieldView(DMField field,PetscViewer viewer)
 72: {
 73:   PetscErrorCode    ierr;
 74:   PetscBool         iascii;

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

 98: /*@C
 99:    DMFieldSetType - set the DMField implementation

101:    Collective on field

103:    Input Parameters:
104: +  field - the DMField context
105: -  type - a known method

107:    Notes:
108:    See "include/petscvec.h" for available methods (for instance)
109: +    DMFIELDDA    - a field defined only by its values at the corners of a DMDA
110: .    DMFIELDDS    - a field defined by a discretization over a mesh set with DMSetField()
111: -    DMFIELDSHELL - a field defined by arbitrary callbacks

113:   Level: advanced

115: .seealso: DMFieldType,
116: @*/
117: PetscErrorCode DMFieldSetType(DMField field,DMFieldType type)
118: {
119:   PetscErrorCode ierr,(*r)(DMField);
120:   PetscBool      match;


126:   PetscObjectTypeCompare((PetscObject)field,type,&match);
127:   if (match) return(0);

129:   PetscFunctionListFind(DMFieldList,type,&r);
130:   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested DMField type %s",type);
131:   /* Destroy the previous private DMField context */
132:   if (field->ops->destroy) {
133:     (*(field)->ops->destroy)(field);
134:   }
135:   PetscMemzero(field->ops,sizeof(*field->ops));
136:   PetscObjectChangeTypeName((PetscObject)field,type);
137:   field->ops->create = r;
138:   (*r)(field);
139:   return(0);
140: }

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

145:   Not Collective

147:   Input Parameter:
148: . field  - The DMField context

150:   Output Parameter:
151: . type - The DMField type name

153:   Level: advanced

155: .seealso: DMFieldSetType()
156: @*/
157: PetscErrorCode  DMFieldGetType(DMField field, DMFieldType *type)
158: {

164:   DMFieldRegisterAll();
165:   *type = ((PetscObject)field)->type_name;
166:   return(0);
167: }

169: /*@
170:   DMFieldGetNumComponents - Returns the number of components in the field

172:   Not collective

174:   Input Parameter:
175: . field - The DMField object

177:   Output Parameter:
178: . nc - The number of field components

180:   Level: intermediate

182: .seealso: DMFieldEvaluate()
183: @*/
184: PetscErrorCode DMFieldGetNumComponents(DMField field, PetscInt *nc)
185: {
189:   *nc = field->numComponents;
190:   return(0);
191: }

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

196:   Not collective

198:   Input Parameter:
199: . field - The DMField object

201:   Output Parameter:
202: . dm - The DM object

204:   Level: intermediate

206: .seealso: DMFieldEvaluate()
207: @*/
208: PetscErrorCode DMFieldGetDM(DMField field, DM *dm)
209: {
213:   *dm = field->dm;
214:   return(0);
215: }

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

220:   Collective on points

222:   Input Parameter:
223: + field - The DMField object
224: . points - The points at which to evaluate the field.  Should have size d x n,
225:            where d is the coordinate dimension of the manifold and n is the number
226:            of points
227: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
228:              If the field is complex and datatype is PETSC_REAL, the real part of the
229:              field is returned.


232:   Output Parameter:
233: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
234:       If B is not NULL, the values of the field are written in this array, varying first by component,
235:       then by point.
236: . D - pointer to data of size d * c * n * sizeof(datatype).
237:       If D is not NULL, the values of the field's spatial derivatives are written in this array,
238:       varying first by the partial derivative component, then by field component, then by point.
239: - H - pointer to data of size d * d * c * n * sizeof(datatype).
240:       If H is not NULL, the values of the field's second spatial derivatives are written in this array,
241:       varying first by the second partial derivative component, then by field component, then by point.

243:   Level: intermediate

245: .seealso: DMFieldGetDM(), DMFieldGetNumComponents(), DMFieldEvaluateFE(), DMFieldEvaluateFV()
246: @*/
247: PetscErrorCode DMFieldEvaluate(DMField field, Vec points, PetscDataType datatype, void *B, void *D, void *H)
248: {

257:   if (field->ops->evaluate) {
258:     (*field->ops->evaluate) (field, points, datatype, B, D, H);
259:   } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
260:   return(0);
261: }

263: /*@
264:   DMFieldEvaluateFE - Evaluate the field and its derivatives on a set of points mapped from
265:   quadrature points on a reference point.  The derivatives are taken with respect to the
266:   reference coordinates.

268:   Not collective

270:   Input Parameter:
271: + field - The DMField object
272: . cellIS - Index set for cells on which to evaluate the field
273: . points - The quadature containing the points in the reference cell at which to evaluate the field.
274: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
275:              If the field is complex and datatype is PETSC_REAL, the real part of the
276:              field is returned.


279:   Output Parameter:
280: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
281:       If B is not NULL, the values of the field are written in this array, varying first by component,
282:       then by point.
283: . D - pointer to data of size d * c * n * sizeof(datatype).
284:       If D is not NULL, the values of the field's spatial derivatives are written in this array,
285:       varying first by the partial derivative component, then by field component, then by point.
286: - H - pointer to data of size d * d * c * n * sizeof(datatype).
287:       If H is not NULL, the values of the field's second spatial derivatives are written in this array,
288:       varying first by the second partial derivative component, then by field component, then by point.

290:   Level: intermediate

292: .seealso: DMFieldGetNumComponents(), DMFieldEvaluate(), DMFieldEvaluateFV()
293: @*/
294: PetscErrorCode DMFieldEvaluateFE(DMField field, IS cellIS, PetscQuadrature points, PetscDataType datatype, void *B, void *D, void *H)
295: {

305:   if (field->ops->evaluateFE) {
306:     (*field->ops->evaluateFE) (field, cellIS, points, datatype, B, D, H);
307:   } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
308:   return(0);
309: }

311: /*@
312:   DMFieldEvaluateFV - Evaluate the mean of a field and its finite volume derivatives on a set of points.

314:   Not collective

316:   Input Parameter:
317: + field - The DMField object
318: . cellIS - Index set for cells on which to evaluate the field
319: - datatype - The PetscDataType of the output arrays: either PETSC_REAL or PETSC_SCALAR.
320:              If the field is complex and datatype is PETSC_REAL, the real part of the
321:              field is returned.


324:   Output Parameter:
325: + B - pointer to data of size c * n * sizeof(datatype), where c is the number of components in the field.
326:       If B is not NULL, the values of the field are written in this array, varying first by component,
327:       then by point.
328: . D - pointer to data of size d * c * n * sizeof(datatype).
329:       If D is not NULL, the values of the field's spatial derivatives are written in this array,
330:       varying first by the partial derivative component, then by field component, then by point.
331: - H - pointer to data of size d * d * c * n * sizeof(datatype).
332:       If H is not NULL, the values of the field's second spatial derivatives are written in this array,
333:       varying first by the second partial derivative component, then by field component, then by point.

335:   Level: intermediate

337: .seealso: DMFieldGetNumComponents(), DMFieldEvaluate(), DMFieldEvaluateFE()
338: @*/
339: PetscErrorCode DMFieldEvaluateFV(DMField field, IS cellIS, PetscDataType datatype, void *B, void *D, void *H)
340: {

349:   if (field->ops->evaluateFV) {
350:     (*field->ops->evaluateFV) (field, cellIS, datatype, B, D, H);
351:   } else SETERRQ (PetscObjectComm((PetscObject)field),PETSC_ERR_SUP,"Not implemented for this type");
352:   return(0);
353: }

355: /*@
356:   DMFieldGetDegree - Get the polynomial degree of a field when pulled back onto the
357:   reference element

359:   Not collective

361:   Input Arguments:
362: + field - the DMField object
363: - cellIS - the index set of points over which we want know the invariance

365:   Output Arguments:
366: + minDegree - the degree of the largest polynomial space contained in the field on each element
367: - maxDegree - the largest degree of the smallest polynomial space containing the field on any element

369:   Level: intermediate

371: .seealso: DMFieldEvaluateFE()
372: @*/
373: PetscErrorCode DMFieldGetDegree(DMField field, IS cellIS, PetscInt *minDegree, PetscInt *maxDegree)
374: {


383:   if (minDegree) *minDegree = -1;
384:   if (maxDegree) *maxDegree = PETSC_MAX_INT;

386:   if (field->ops->getDegree) {
387:     (*field->ops->getDegree) (field,cellIS,minDegree,maxDegree);
388:   }
389:   return(0);
390: }

392: /*@
393:   DMFieldCreateDefaultQuadrature - Creates a quadrature sufficient to integrate the field on the selected
394:   points via pullback onto the reference element

396:   Not collective

398:   Input Arguments:
399: + field - the DMField object
400: - pointIS - the index set of points over which we wish to integrate the field

402:   Output Arguments:
403: . quad - a PetscQuadrature object

405:   Level: developer

407: .seealso: DMFieldEvaluteFE(), DMFieldGetDegree()
408: @*/
409: PetscErrorCode DMFieldCreateDefaultQuadrature(DMField field, IS pointIS, PetscQuadrature *quad)
410: {


418:   *quad = NULL;
419:   if (field->ops->createDefaultQuadrature) {
420:     (*field->ops->createDefaultQuadrature)(field, pointIS, quad);
421:   }
422:   return(0);
423: }

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

428:   Not collective

430:   Input Arguments:
431: + field - the DMField object
432: . pointIS - the index set of points over which we wish to integrate the field
433: . quad - the quadrature points at which to evaluate the geometric factors
434: - faceData - whether additional data for facets (the normal vectors and adjacent cells) should
435:   be calculated

437:   Output Arguments:
438: . geom - the geometric factors

440:   Level: developer

442: .seealso: DMFieldEvaluateFE(), DMFieldCreateDefaulteQuadrature(), DMFieldGetDegree()
443: @*/
444: PetscErrorCode DMFieldCreateFEGeom(DMField field, IS pointIS, PetscQuadrature quad, PetscBool faceData, PetscFEGeom **geom)
445: {
446:   PetscInt       dim, dE;
447:   PetscInt       nPoints;
448:   PetscInt       maxDegree;
449:   PetscFEGeom    *g;

456:   ISGetLocalSize(pointIS,&nPoints);
457:   dE = field->numComponents;
458:   PetscFEGeomCreate(quad,nPoints,dE,faceData,&g);
459:   DMFieldEvaluateFE(field,pointIS,quad,PETSC_REAL,g->v,g->J,NULL);
460:   dim = g->dim;
461:   if (dE > dim) {
462:     /* space out J and make square Jacobians */
463:     PetscInt  i, j, k, N = g->numPoints * g->numCells;

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

468:       for (j = 0; j < dE; j++) {
469:         for (k = 0; k < dim; k++) {
470:           J[j*dE + k] = g->J[i*dE*dim + j*dim + k];
471:         }
472:       }
473:       switch (dim) {
474:       case 0:
475:         for (j = 0; j < dE; j++) {
476:           for (k = 0; k < dE; k++) {
477:             J[j * dE + k] = (j == k) ? 1. : 0.;
478:           }
479:         }
480:         break;
481:       case 1:
482:         if (dE == 2) {
483:           PetscReal norm = PetscSqrtReal(J[0] * J[0] + J[2] * J[2]);

485:           J[1] = -J[2] / norm;
486:           J[3] =  J[0] / norm;
487:         } else {
488:           PetscReal inorm = 1./PetscSqrtReal(J[0] * J[0] + J[3] * J[3] + J[6] * J[6]);
489:           PetscReal x = J[0] * inorm;
490:           PetscReal y = J[3] * inorm;
491:           PetscReal z = J[6] * inorm;

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

496:             J[1] = -y;              J[2] = -z;
497:             J[4] = 1. - y*y*inv1pX; J[5] =     -y*z*inv1pX;
498:             J[7] =     -y*z*inv1pX; J[8] = 1. - z*z*inv1pX;
499:           } else {
500:             PetscReal inv1mX   = 1./ (1. - x);

502:             J[1] = z;               J[2] = y;
503:             J[4] =     -y*z*inv1mX; J[5] = 1. - y*y*inv1mX;
504:             J[7] = 1. - z*z*inv1mX; J[8] =     -y*z*inv1mX;
505:           }
506:         }
507:         break;
508:       case 2:
509:         {
510:           PetscReal inorm;

512:           J[2] = J[3] * J[7] - J[6] * J[4];
513:           J[5] = J[6] * J[1] - J[0] * J[7];
514:           J[8] = J[0] * J[4] - J[3] * J[1];

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

518:           J[2] *= inorm;
519:           J[5] *= inorm;
520:           J[8] *= inorm;
521:         }
522:         break;
523:       }
524:       for (j = 0; j < dE*dE; j++) {
525:         g->J[i*dE*dE + j] = J[j];
526:       }
527:     }
528:   }
529:   PetscFEGeomComplete(g);
530:   DMFieldGetDegree(field,pointIS,NULL,&maxDegree);
531:   g->isAffine = (maxDegree <= 1) ? PETSC_TRUE : PETSC_FALSE;
532:   if (faceData) {
533:     if (!field->ops->computeFaceData) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "DMField implementation does not compute face data\n");
534:     (*field->ops->computeFaceData) (field, pointIS, quad, g);
535:   }
536:   *geom = g;
537:   return(0);
538: }