1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
2 #include <petsc/private/dtimpl.h> /*I "petscdt.h" I*/
3
PetscSpacePointView_Ascii(PetscSpace sp,PetscViewer viewer)4 static PetscErrorCode PetscSpacePointView_Ascii(PetscSpace sp, PetscViewer viewer)
5 {
6 PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;
7 PetscViewerFormat format;
8
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 }
19
PetscSpaceView_Point(PetscSpace sp,PetscViewer viewer)20 static PetscErrorCode PetscSpaceView_Point(PetscSpace sp, PetscViewer viewer)
21 {
22 PetscBool isascii;
23
24 PetscFunctionBegin;
25 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
26 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
27 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
28 if (isascii) PetscCall(PetscSpacePointView_Ascii(sp, viewer));
29 PetscFunctionReturn(PETSC_SUCCESS);
30 }
31
PetscSpaceSetUp_Point(PetscSpace sp)32 static PetscErrorCode PetscSpaceSetUp_Point(PetscSpace sp)
33 {
34 PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;
35
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 }
43
PetscSpaceDestroy_Point(PetscSpace sp)44 static PetscErrorCode PetscSpaceDestroy_Point(PetscSpace sp)
45 {
46 PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;
47
48 PetscFunctionBegin;
49 PetscCall(PetscQuadratureDestroy(&pt->quad));
50 PetscCall(PetscFree(pt));
51 PetscFunctionReturn(PETSC_SUCCESS);
52 }
53
PetscSpaceGetDimension_Point(PetscSpace sp,PetscInt * dim)54 static PetscErrorCode PetscSpaceGetDimension_Point(PetscSpace sp, PetscInt *dim)
55 {
56 PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;
57
58 PetscFunctionBegin;
59 *dim = pt->quad->numPoints;
60 PetscFunctionReturn(PETSC_SUCCESS);
61 }
62
PetscSpaceEvaluate_Point(PetscSpace sp,PetscInt npoints,const PetscReal points[],PetscReal B[],PetscReal D[],PetscReal H[])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;
67
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 }
92
PetscSpaceInitialize_Point(PetscSpace sp)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 }
104
105 /*MC
106 PETSCSPACEPOINT = "point" - A `PetscSpace` object that encapsulates functions defined on a set of quadrature points.
107
108 Level: intermediate
109
110 .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
111 M*/
112
PetscSpaceCreate_Point(PetscSpace sp)113 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Point(PetscSpace sp)
114 {
115 PetscSpace_Point *pt;
116
117 PetscFunctionBegin;
118 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
119 PetscCall(PetscNew(&pt));
120 sp->data = pt;
121
122 sp->Nv = 0;
123 sp->maxDegree = PETSC_INT_MAX;
124 PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, &pt->quad));
125 PetscCall(PetscQuadratureSetData(pt->quad, 0, 1, 0, NULL, NULL));
126
127 PetscCall(PetscSpaceInitialize_Point(sp));
128 PetscFunctionReturn(PETSC_SUCCESS);
129 }
130
131 /*@
132 PetscSpacePointSetPoints - Sets the evaluation points for the space to coincide with the points of a quadrature rule
133
134 Logically Collective
135
136 Input Parameters:
137 + sp - The `PetscSpace`
138 - q - The `PetscQuadrature` defining the points
139
140 Level: intermediate
141
142 .seealso: `PetscSpace`, `PetscQuadrature`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
143 @*/
PetscSpacePointSetPoints(PetscSpace sp,PetscQuadrature q)144 PetscErrorCode PetscSpacePointSetPoints(PetscSpace sp, PetscQuadrature q)
145 {
146 PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;
147
148 PetscFunctionBegin;
149 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
150 PetscValidHeaderSpecific(q, PETSCQUADRATURE_CLASSID, 2);
151 PetscCall(PetscQuadratureDestroy(&pt->quad));
152 PetscCall(PetscQuadratureDuplicate(q, &pt->quad));
153 PetscFunctionReturn(PETSC_SUCCESS);
154 }
155
156 /*@
157 PetscSpacePointGetPoints - Gets the evaluation points for the space as the points of a quadrature rule
158
159 Logically Collective
160
161 Input Parameter:
162 . sp - The `PetscSpace`
163
164 Output Parameter:
165 . q - The `PetscQuadrature` defining the points
166
167 Level: intermediate
168
169 .seealso: `PetscSpace`, `PetscQuadrature`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
170 @*/
PetscSpacePointGetPoints(PetscSpace sp,PetscQuadrature * q)171 PetscErrorCode PetscSpacePointGetPoints(PetscSpace sp, PetscQuadrature *q)
172 {
173 PetscSpace_Point *pt = (PetscSpace_Point *)sp->data;
174
175 PetscFunctionBegin;
176 PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
177 PetscAssertPointer(q, 2);
178 *q = pt->quad;
179 PetscFunctionReturn(PETSC_SUCCESS);
180 }
181