Actual source code: dmfield.c

petsc-3.11.2 2019-05-18
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 DMField

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: .keywords: DMField, set, type

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


128:   PetscObjectTypeCompare((PetscObject)field,type,&match);
129:   if (match) return(0);

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

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

147:   Not Collective

149:   Input Parameter:
150: . field  - The DMField context

152:   Output Parameter:
153: . type - The DMField type name

155:   Level: advanced

157: .keywords: DMField, get, type, name
158: .seealso: DMFieldSetType()
159: @*/
160: PetscErrorCode  DMFieldGetType(DMField field, DMFieldType *type)
161: {

167:   DMFieldRegisterAll();
168:   *type = ((PetscObject)field)->type_name;
169:   return(0);
170: }

172: /*@
173:   DMFieldGetNumComponents - Returns the number of components in the field

175:   Not collective

177:   Input Parameter:
178: . field - The DMField object

180:   Output Parameter:
181: . nc - The number of field components

183:   Level: intermediate

185: .seealso: DMFieldEvaluate()
186: @*/
187: PetscErrorCode DMFieldGetNumComponents(DMField field, PetscInt *nc)
188: {
192:   *nc = field->numComponents;
193:   return(0);
194: }

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

199:   Not collective

201:   Input Parameter:
202: . field - The DMField object

204:   Output Parameter:
205: . dm - The DM object

207:   Level: intermediate

209: .seealso: DMFieldEvaluate()
210: @*/
211: PetscErrorCode DMFieldGetDM(DMField field, DM *dm)
212: {
216:   *dm = field->dm;
217:   return(0);
218: }

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

223:   Collective on Vec

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


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

246:   Level: intermediate

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

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

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

271:   Not collective

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


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

293:   Level: intermediate

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

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

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

317:   Not collective

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


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

338:   Level: intermediate

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

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

358: /*@
359:   DMFieldGetDegree - Get the polynomial degree of a field when pulled back onto the
360:   reference element

362:   Not collective

364:   Input Arguments:
365: + field - the DMField object
366: - cellIS - the index set of points over which we want know the invariance

368:   Output Arguments:
369: + minDegree - the degree of the largest polynomial space contained in the field on each element
370: - maxDegree - the largest degree of the smallest polynomial space containing the field on any element

372:   Level: intermediate

374: .seealso: DMFieldEvaluateFE()
375: @*/
376: PetscErrorCode DMFieldGetDegree(DMField field, IS cellIS, PetscInt *minDegree, PetscInt *maxDegree)
377: {


386:   if (minDegree) *minDegree = -1;
387:   if (maxDegree) *maxDegree = PETSC_MAX_INT;

389:   if (field->ops->getDegree) {
390:     (*field->ops->getDegree) (field,cellIS,minDegree,maxDegree);
391:   }
392:   return(0);
393: }

395: /*@
396:   DMFieldCreateDefaultQuadrature - Creates a quadrature sufficient to integrate the field on the selected
397:   points via pullback onto the reference element

399:   Not collective

401:   Input Arguments:
402: + field - the DMField object
403: - pointIS - the index set of points over which we wish to integrate the field

405:   Output Arguments:
406: . quad - a PetscQuadrature object

408:   Level: developer

410: .seealso: DMFieldEvaluteFE(), DMFieldGetDegree()
411: @*/
412: PetscErrorCode DMFieldCreateDefaultQuadrature(DMField field, IS pointIS, PetscQuadrature *quad)
413: {


421:   *quad = NULL;
422:   if (field->ops->createDefaultQuadrature) {
423:     (*field->ops->createDefaultQuadrature)(field, pointIS, quad);
424:   }
425:   return(0);
426: }

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

431:   Not collective

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

440:   Output Arguments:
441: . geom - the geometric factors

443:   Level: developer

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

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

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

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

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

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

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

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

515:           J[2] = J[3] * J[7] - J[6] * J[4];
516:           J[5] = J[6] * J[1] - J[0] * J[7];
517:           J[8] = J[0] * J[4] - J[3] * J[1];

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

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