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:   PetscFunctionBegin;
 10:   PetscCall(PetscViewerGetFormat(viewer, &format));
 11:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
 12:     PetscCall(PetscViewerASCIIPrintf(viewer, "Point space in dimension %" PetscInt_FMT ":\n", sp->Nv));
 13:     PetscCall(PetscViewerASCIIPushTab(viewer));
 14:     PetscCall(PetscQuadratureView(pt->quad, viewer));
 15:     PetscCall(PetscViewerASCIIPopTab(viewer));
 16:   } else PetscCall(PetscViewerASCIIPrintf(viewer, "Point space in dimension %" PetscInt_FMT " on %" PetscInt_FMT " points\n", sp->Nv, pt->quad->numPoints));
 17:   PetscFunctionReturn(PETSC_SUCCESS);
 18: }

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

 24:   PetscFunctionBegin;
 27:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
 28:   if (iascii) PetscCall(PetscSpacePointView_Ascii(sp, viewer));
 29:   PetscFunctionReturn(PETSC_SUCCESS);
 30: }

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

 36:   PetscFunctionBegin;
 37:   if (!pt->quad->points && sp->degree >= 0) {
 38:     PetscCall(PetscQuadratureDestroy(&pt->quad));
 39:     PetscCall(PetscDTStroudConicalQuadrature(sp->Nv, sp->Nc, PetscMax(sp->degree + 1, 1), -1.0, 1.0, &pt->quad));
 40:   }
 41:   PetscFunctionReturn(PETSC_SUCCESS);
 42: }

 44: static PetscErrorCode PetscSpaceDestroy_Point(PetscSpace sp)
 45: {
 46:   PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;

 48:   PetscFunctionBegin;
 49:   PetscCall(PetscQuadratureDestroy(&pt->quad));
 50:   PetscCall(PetscFree(pt));
 51:   PetscFunctionReturn(PETSC_SUCCESS);
 52: }

 54: static PetscErrorCode PetscSpaceGetDimension_Point(PetscSpace sp, PetscInt *dim)
 55: {
 56:   PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;

 58:   PetscFunctionBegin;
 59:   *dim = pt->quad->numPoints;
 60:   PetscFunctionReturn(PETSC_SUCCESS);
 61: }

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

 68:   PetscFunctionBegin;
 69:   PetscCheck(npoints == pt->quad->numPoints, PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot evaluate Point space on %" PetscInt_FMT " points != %" PetscInt_FMT " size", npoints, pt->quad->numPoints);
 70:   PetscCall(PetscArrayzero(B, npoints * pdim));
 71:   for (p = 0; p < npoints; ++p) {
 72:     for (i = 0; i < pdim; ++i) {
 73:       for (d = 0; d < dim; ++d) {
 74:         if (PetscAbsReal(points[p * dim + d] - pt->quad->points[p * dim + d]) > 1.0e-10) break;
 75:       }
 76:       if (d >= dim) {
 77:         B[p * pdim + i] = 1.0;
 78:         break;
 79:       }
 80:     }
 81:   }
 82:   /* Replicate for other components */
 83:   for (c = 1; c < sp->Nc; ++c) {
 84:     for (p = 0; p < npoints; ++p) {
 85:       for (i = 0; i < pdim; ++i) B[(c * npoints + p) * pdim + i] = B[p * pdim + i];
 86:     }
 87:   }
 88:   if (D) PetscCall(PetscArrayzero(D, npoints * pdim * dim));
 89:   if (H) PetscCall(PetscArrayzero(H, npoints * pdim * dim * dim));
 90:   PetscFunctionReturn(PETSC_SUCCESS);
 91: }

 93: static PetscErrorCode PetscSpaceInitialize_Point(PetscSpace sp)
 94: {
 95:   PetscFunctionBegin;
 96:   sp->ops->setfromoptions = NULL;
 97:   sp->ops->setup          = PetscSpaceSetUp_Point;
 98:   sp->ops->view           = PetscSpaceView_Point;
 99:   sp->ops->destroy        = PetscSpaceDestroy_Point;
100:   sp->ops->getdimension   = PetscSpaceGetDimension_Point;
101:   sp->ops->evaluate       = PetscSpaceEvaluate_Point;
102:   PetscFunctionReturn(PETSC_SUCCESS);
103: }

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

108:   Level: intermediate

110: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
111: M*/

113: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Point(PetscSpace sp)
114: {
115:   PetscSpace_Point *pt;

117:   PetscFunctionBegin;
119:   PetscCall(PetscNew(&pt));
120:   sp->data = pt;

122:   sp->Nv        = 0;
123:   sp->maxDegree = PETSC_MAX_INT;
124:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &pt->quad));
125:   PetscCall(PetscQuadratureSetData(pt->quad, 0, 1, 0, NULL, NULL));

127:   PetscCall(PetscSpaceInitialize_Point(sp));
128:   PetscFunctionReturn(PETSC_SUCCESS);
129: }

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

134:   Logically Collective

136:   Input Parameters:
137: + sp - The `PetscSpace`
138: - q  - The `PetscQuadrature` defining the points

140:   Level: intermediate

142: .seealso: `PetscSpace`, `PetscQuadrature`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
143: @*/
144: PetscErrorCode PetscSpacePointSetPoints(PetscSpace sp, PetscQuadrature q)
145: {
146:   PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;

148:   PetscFunctionBegin;
151:   PetscCall(PetscQuadratureDestroy(&pt->quad));
152:   PetscCall(PetscQuadratureDuplicate(q, &pt->quad));
153:   PetscFunctionReturn(PETSC_SUCCESS);
154: }

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

159:   Logically Collective

161:   Input Parameter:
162: . sp - The `PetscSpace`

164:   Output Parameter:
165: . q - The `PetscQuadrature` defining the points

167:   Level: intermediate

169: .seealso: `PetscSpace`, `PetscQuadrature`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
170: @*/
171: PetscErrorCode PetscSpacePointGetPoints(PetscSpace sp, PetscQuadrature *q)
172: {
173:   PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;

175:   PetscFunctionBegin;
177:   PetscAssertPointer(q, 2);
178:   *q = pt->quad;
179:   PetscFunctionReturn(PETSC_SUCCESS);
180: }