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 %d:\n", sp->Nv);
 12:     PetscViewerASCIIPushTab(viewer);
 13:     PetscQuadratureView(pt->quad, viewer);
 14:     PetscViewerASCIIPopTab(viewer);
 15:   } else {
 16:     PetscViewerASCIIPrintf(viewer, "Point space in dimension %d on %d points\n", sp->Nv, pt->quad->numPoints);
 17:   }
 18:   return 0;
 19: }

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

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

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

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

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

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

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

 56:   *dim = pt->quad->numPoints;
 57:   return 0;
 58: }

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

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

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

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

102:   Level: intermediate

104: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
105: M*/

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

112:   PetscNewLog(sp,&pt);
113:   sp->data = pt;

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

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

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

127:   Logically collective

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

133:   Level: intermediate

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

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

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

151:   Logically collective

153:   Input Parameter:
154: . sp - The PetscSpace

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

159:   Level: intermediate

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

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