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: }