Actual source code: spacepoint.c
petsc-3.10.5 2019-03-28
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: }