Actual source code: spacepoint.c

  1: #include <petsc/private/petscfeimpl.h>
  2: #include <petsc/private/dtimpl.h>

  4: static PetscErrorCode PetscSpacePointView_Ascii(PetscSpace sp, PetscViewer viewer)
  5: {
  6:   PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;
  7:   PetscViewerFormat format;

  9:   PetscViewerGetFormat(viewer, &format);
 10:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
 11:     PetscViewerASCIIPrintf(viewer, "Point space in dimension %" PetscInt_FMT ":\n", sp->Nv);
 12:     PetscViewerASCIIPushTab(viewer);
 13:     PetscQuadratureView(pt->quad, viewer);
 14:     PetscViewerASCIIPopTab(viewer);
 15:   } else PetscViewerASCIIPrintf(viewer, "Point space in dimension %" PetscInt_FMT " on %" PetscInt_FMT " points\n", sp->Nv, pt->quad->numPoints);
 16:   return 0;
 17: }

 19: static PetscErrorCode PetscSpaceView_Point(PetscSpace sp, PetscViewer viewer)
 20: {
 21:   PetscBool iascii;

 25:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
 26:   if (iascii) PetscSpacePointView_Ascii(sp, viewer);
 27:   return 0;
 28: }

 30: static PetscErrorCode PetscSpaceSetUp_Point(PetscSpace sp)
 31: {
 32:   PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;

 34:   if (!pt->quad->points && sp->degree >= 0) {
 35:     PetscQuadratureDestroy(&pt->quad);
 36:     PetscDTStroudConicalQuadrature(sp->Nv, sp->Nc, PetscMax(sp->degree + 1, 1), -1.0, 1.0, &pt->quad);
 37:   }
 38:   return 0;
 39: }

 41: static PetscErrorCode PetscSpaceDestroy_Point(PetscSpace sp)
 42: {
 43:   PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;

 45:   PetscQuadratureDestroy(&pt->quad);
 46:   PetscFree(pt);
 47:   return 0;
 48: }

 50: static PetscErrorCode PetscSpaceGetDimension_Point(PetscSpace sp, PetscInt *dim)
 51: {
 52:   PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;

 54:   *dim = pt->quad->numPoints;
 55:   return 0;
 56: }

 58: static PetscErrorCode PetscSpaceEvaluate_Point(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
 59: {
 60:   PetscSpace_Point *pt  = (PetscSpace_Point *)sp->data;
 61:   PetscInt          dim = sp->Nv, pdim = pt->quad->numPoints, d, p, i, c;

 64:   PetscArrayzero(B, npoints * pdim);
 65:   for (p = 0; p < npoints; ++p) {
 66:     for (i = 0; i < pdim; ++i) {
 67:       for (d = 0; d < dim; ++d) {
 68:         if (PetscAbsReal(points[p * dim + d] - pt->quad->points[p * dim + d]) > 1.0e-10) break;
 69:       }
 70:       if (d >= dim) {
 71:         B[p * pdim + i] = 1.0;
 72:         break;
 73:       }
 74:     }
 75:   }
 76:   /* Replicate for other components */
 77:   for (c = 1; c < sp->Nc; ++c) {
 78:     for (p = 0; p < npoints; ++p) {
 79:       for (i = 0; i < pdim; ++i) B[(c * npoints + p) * pdim + i] = B[p * pdim + i];
 80:     }
 81:   }
 82:   if (D) PetscArrayzero(D, npoints * pdim * dim);
 83:   if (H) PetscArrayzero(H, npoints * pdim * dim * dim);
 84:   return 0;
 85: }

 87: static PetscErrorCode PetscSpaceInitialize_Point(PetscSpace sp)
 88: {
 89:   sp->ops->setfromoptions = NULL;
 90:   sp->ops->setup          = PetscSpaceSetUp_Point;
 91:   sp->ops->view           = PetscSpaceView_Point;
 92:   sp->ops->destroy        = PetscSpaceDestroy_Point;
 93:   sp->ops->getdimension   = PetscSpaceGetDimension_Point;
 94:   sp->ops->evaluate       = PetscSpaceEvaluate_Point;
 95:   return 0;
 96: }

 98: /*MC
 99:   PETSCSPACEPOINT = "point" - A `PetscSpace` object that encapsulates functions defined on a set of quadrature points.

101:   Level: intermediate

103: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
104: M*/

106: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Point(PetscSpace sp)
107: {
108:   PetscSpace_Point *pt;

111:   PetscNew(&pt);
112:   sp->data = pt;

114:   sp->Nv        = 0;
115:   sp->maxDegree = PETSC_MAX_INT;
116:   PetscQuadratureCreate(PETSC_COMM_SELF, &pt->quad);
117:   PetscQuadratureSetData(pt->quad, 0, 1, 0, NULL, NULL);

119:   PetscSpaceInitialize_Point(sp);
120:   return 0;
121: }

123: /*@
124:   PetscSpacePointSetPoints - Sets the evaluation points for the space to coincide with the points of a quadrature rule

126:   Logically collective

128:   Input Parameters:
129: + sp - The `PetscSpace`
130: - q  - The `PetscQuadrature` defining the points

132:   Level: intermediate

134: .seealso: `PetscSpace`, `PetscQuadrature`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
135: @*/
136: PetscErrorCode PetscSpacePointSetPoints(PetscSpace sp, PetscQuadrature q)
137: {
138:   PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;

142:   PetscQuadratureDestroy(&pt->quad);
143:   PetscQuadratureDuplicate(q, &pt->quad);
144:   return 0;
145: }

147: /*@
148:   PetscSpacePointGetPoints - Gets the evaluation points for the space as the points of a quadrature rule

150:   Logically collective

152:   Input Parameter:
153: . sp - The `PetscSpace`

155:   Output Parameter:
156: . q  - The `PetscQuadrature` defining the points

158:   Level: intermediate

160: .seealso: `PetscSpace`, `PetscQuadrature`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
161: @*/
162: PetscErrorCode PetscSpacePointGetPoints(PetscSpace sp, PetscQuadrature *q)
163: {
164:   PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;

168:   *q = pt->quad;
169:   return 0;
170: }