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

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

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

 23: PetscErrorCode PetscSpaceView_Point(PetscSpace sp, PetscViewer viewer)
 24: {
 25:   PetscBool      iascii;

 31:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
 32:   if (iascii) {PetscSpacePointView_Ascii(sp, viewer);}
 33:   return(0);
 34: }

 36: PetscErrorCode PetscSpaceSetUp_Point(PetscSpace sp)
 37: {
 38:   PetscSpace_Point *pt = (PetscSpace_Point *) sp->data;
 39:   PetscErrorCode    ierr;

 42:   if (!pt->quad->points && sp->degree >= 0) {
 43:     PetscQuadratureDestroy(&pt->quad);
 44:     PetscDTGaussJacobiQuadrature(sp->Nv, sp->Nc, PetscMax(sp->degree + 1, 1), -1.0, 1.0, &pt->quad);
 45:   }
 46:   return(0);
 47: }

 49: PetscErrorCode PetscSpaceDestroy_Point(PetscSpace sp)
 50: {
 51:   PetscSpace_Point *pt = (PetscSpace_Point *) sp->data;
 52:   PetscErrorCode    ierr;

 55:   PetscQuadratureDestroy(&pt->quad);
 56:   PetscFree(pt);
 57:   return(0);
 58: }

 60: PetscErrorCode PetscSpaceGetDimension_Point(PetscSpace sp, PetscInt *dim)
 61: {
 62:   PetscSpace_Point *pt = (PetscSpace_Point *) sp->data;

 65:   *dim = pt->quad->numPoints;
 66:   return(0);
 67: }

 69: PetscErrorCode PetscSpaceEvaluate_Point(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
 70: {
 71:   PetscSpace_Point *pt  = (PetscSpace_Point *) sp->data;
 72:   PetscInt          dim = sp->Nv, pdim = pt->quad->numPoints, d, p, i, c;
 73:   PetscErrorCode    ierr;

 76:   if (npoints != pt->quad->numPoints) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot evaluate Point space on %d points != %d size", npoints, pt->quad->numPoints);
 77:   PetscMemzero(B, npoints*pdim * sizeof(PetscReal));
 78:   for (p = 0; p < npoints; ++p) {
 79:     for (i = 0; i < pdim; ++i) {
 80:       for (d = 0; d < dim; ++d) {
 81:         if (PetscAbsReal(points[p*dim+d] - pt->quad->points[p*dim+d]) > 1.0e-10) break;
 82:       }
 83:       if (d >= dim) {B[p*pdim+i] = 1.0; break;}
 84:     }
 85:   }
 86:   /* Replicate for other components */
 87:   for (c = 1; c < sp->Nc; ++c) {
 88:     for (p = 0; p < npoints; ++p) {
 89:       for (i = 0; i < pdim; ++i) {
 90:         B[(c*npoints + p)*pdim + i] = B[p*pdim + i];
 91:       }
 92:     }
 93:   }
 94:   if (D) {PetscMemzero(D, npoints*pdim*dim * sizeof(PetscReal));}
 95:   if (H) {PetscMemzero(H, npoints*pdim*dim*dim * sizeof(PetscReal));}
 96:   return(0);
 97: }

 99: PetscErrorCode PetscSpaceInitialize_Point(PetscSpace sp)
100: {
102:   sp->ops->setfromoptions = NULL;
103:   sp->ops->setup          = PetscSpaceSetUp_Point;
104:   sp->ops->view           = PetscSpaceView_Point;
105:   sp->ops->destroy        = PetscSpaceDestroy_Point;
106:   sp->ops->getdimension   = PetscSpaceGetDimension_Point;
107:   sp->ops->evaluate       = PetscSpaceEvaluate_Point;
108:   return(0);
109: }

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

114:   Level: intermediate

116: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
117: M*/

119: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Point(PetscSpace sp)
120: {
121:   PetscSpace_Point *pt;
122:   PetscErrorCode    ierr;

126:   PetscNewLog(sp,&pt);
127:   sp->data = pt;

129:   sp->Nv = 0;
130:   sp->maxDegree = PETSC_MAX_INT;
131:   PetscQuadratureCreate(PETSC_COMM_SELF, &pt->quad);
132:   PetscQuadratureSetData(pt->quad, 0, 1, 0, NULL, NULL);

134:   PetscSpaceInitialize_Point(sp);
135:   return(0);
136: }

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

141:   Logically collective

143:   Input Parameters:
144: + sp - The PetscSpace
145: - q  - The PetscQuadrature defining the points

147:   Level: intermediate

149: .keywords: PetscSpacePoint
150: .seealso: PetscSpaceCreate(), PetscSpaceSetType()
151: @*/
152: PetscErrorCode PetscSpacePointSetPoints(PetscSpace sp, PetscQuadrature q)
153: {
154:   PetscSpace_Point *pt = (PetscSpace_Point *) sp->data;
155:   PetscErrorCode    ierr;

160:   PetscQuadratureDestroy(&pt->quad);
161:   PetscQuadratureDuplicate(q, &pt->quad);
162:   return(0);
163: }

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

168:   Logically collective

170:   Input Parameter:
171: . sp - The PetscSpace

173:   Output Parameter:
174: . q  - The PetscQuadrature defining the points

176:   Level: intermediate

178: .keywords: PetscSpacePoint
179: .seealso: PetscSpaceCreate(), PetscSpaceSetType()
180: @*/
181: PetscErrorCode PetscSpacePointGetPoints(PetscSpace sp, PetscQuadrature *q)
182: {
183:   PetscSpace_Point *pt = (PetscSpace_Point *) sp->data;

188:   *q = pt->quad;
189:   return(0);
190: }