xref: /petsc/src/dm/dt/space/impls/wxy/spacewxy.c (revision bcda9346efad4e5ba2d553af84eb238771ba1e25) 
1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
2 
PetscSpaceSetFromOptions_WXY(PetscSpace sp,PetscOptionItems PetscOptionsObject)3 static PetscErrorCode PetscSpaceSetFromOptions_WXY(PetscSpace sp, PetscOptionItems PetscOptionsObject)
4 {
5   PetscFunctionBegin;
6   PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace WXY options");
7   PetscOptionsHeadEnd();
8   PetscFunctionReturn(PETSC_SUCCESS);
9 }
10 
PetscSpacePolynomialView_Ascii(PetscSpace sp,PetscViewer v)11 static PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer v)
12 {
13   PetscFunctionBegin;
14   PetscCall(PetscViewerASCIIPrintf(v, "WXY space of degree %" PetscInt_FMT "\n", sp->degree));
15   PetscFunctionReturn(PETSC_SUCCESS);
16 }
17 
PetscSpaceView_WXY(PetscSpace sp,PetscViewer viewer)18 static PetscErrorCode PetscSpaceView_WXY(PetscSpace sp, PetscViewer viewer)
19 {
20   PetscBool isascii;
21 
22   PetscFunctionBegin;
23   PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
24   PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
25   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
26   if (isascii) PetscCall(PetscSpacePolynomialView_Ascii(sp, viewer));
27   PetscFunctionReturn(PETSC_SUCCESS);
28 }
29 
PetscSpaceDestroy_WXY(PetscSpace sp)30 static PetscErrorCode PetscSpaceDestroy_WXY(PetscSpace sp)
31 {
32   PetscSpace_WXY *wxy = (PetscSpace_WXY *)sp->data;
33 
34   PetscFunctionBegin;
35   PetscCall(PetscFree(wxy));
36   PetscFunctionReturn(PETSC_SUCCESS);
37 }
38 
PetscSpaceSetUp_WXY(PetscSpace sp)39 static PetscErrorCode PetscSpaceSetUp_WXY(PetscSpace sp)
40 {
41   PetscSpace_WXY *wxy = (PetscSpace_WXY *)sp->data;
42 
43   PetscFunctionBegin;
44   if (wxy->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
45   PetscCheck(sp->degree >= 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Negative degree %" PetscInt_FMT " invalid", sp->degree);
46   sp->maxDegree    = sp->degree;
47   wxy->setupcalled = PETSC_TRUE;
48   PetscFunctionReturn(PETSC_SUCCESS);
49 }
50 
PetscSpaceGetDimension_WXY(PetscSpace sp,PetscInt * dim)51 static PetscErrorCode PetscSpaceGetDimension_WXY(PetscSpace sp, PetscInt *dim)
52 {
53   PetscFunctionBegin;
54   *dim = 6;
55   PetscFunctionReturn(PETSC_SUCCESS);
56 }
57 
PetscSpaceEvaluate_WXY(PetscSpace sp,PetscInt npoints,const PetscReal points[],PetscReal B[],PetscReal D[],PetscReal H[])58 static PetscErrorCode PetscSpaceEvaluate_WXY(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
59 {
60   PetscSpace_WXY *wxy = (PetscSpace_WXY *)sp->data;
61   PetscInt        dim = sp->Nv;
62   PetscInt        Nb  = 6;
63   PetscInt        Nc  = 3;
64 
65   PetscFunctionBegin;
66   if (!wxy->setupcalled) {
67     PetscCall(PetscSpaceSetUp(sp));
68     PetscCall(PetscSpaceEvaluate(sp, npoints, points, B, D, H));
69     PetscFunctionReturn(PETSC_SUCCESS);
70   }
71   PetscCheck((sp->Nc == 3) && (sp->Nv == 3), PETSC_COMM_SELF, PETSC_ERR_PLIB, "WXY space must have 3 variables and 3 components");
72   if (B) {
73     PetscInt p_inc = Nb * Nc;
74     PetscInt b_inc = Nc;
75     PetscInt c_inc = 1;
76 
77     for (PetscInt p = 0; p < npoints; p++) {
78       const PetscReal x = points[p * dim + 0];
79       const PetscReal y = points[p * dim + 1];
80       const PetscReal z = points[p * dim + 2];
81 
82       /* {2 y z, 0, 0} */
83       B[p * p_inc + 0 * b_inc + 0 * c_inc] = 2. * y * z;
84       B[p * p_inc + 0 * b_inc + 1 * c_inc] = 0.;
85       B[p * p_inc + 0 * b_inc + 2 * c_inc] = 0.;
86       /* {0, 2 x z, 0} */
87       B[p * p_inc + 1 * b_inc + 0 * c_inc] = 0.;
88       B[p * p_inc + 1 * b_inc + 1 * c_inc] = 2. * x * z;
89       B[p * p_inc + 1 * b_inc + 2 * c_inc] = 0.;
90       /* {0, 2 y z, -z^2} */
91       B[p * p_inc + 2 * b_inc + 0 * c_inc] = 0.;
92       B[p * p_inc + 2 * b_inc + 1 * c_inc] = 2. * y * z;
93       B[p * p_inc + 2 * b_inc + 2 * c_inc] = -z * z;
94       /* {2 x z, 0, -z^2} */
95       B[p * p_inc + 3 * b_inc + 0 * c_inc] = 2. * x * z;
96       B[p * p_inc + 3 * b_inc + 1 * c_inc] = 0.;
97       B[p * p_inc + 3 * b_inc + 2 * c_inc] = -z * z;
98       /* {x^2, x y, -3 x z} */
99       B[p * p_inc + 4 * b_inc + 0 * c_inc] = x * x;
100       B[p * p_inc + 4 * b_inc + 1 * c_inc] = x * y;
101       B[p * p_inc + 4 * b_inc + 2 * c_inc] = -3. * x * z;
102       /* {x y, y^2, -3 y z} */
103       B[p * p_inc + 5 * b_inc + 0 * c_inc] = x * y;
104       B[p * p_inc + 5 * b_inc + 1 * c_inc] = y * y;
105       B[p * p_inc + 5 * b_inc + 2 * c_inc] = -3. * y * z;
106     }
107   }
108   if (D) {
109     PetscInt p_inc = Nb * Nc * dim;
110     PetscInt b_inc = Nc * dim;
111     PetscInt c_inc = dim;
112 
113     for (PetscInt p = 0; p < npoints; p++) {
114       const PetscReal x = points[p * dim + 0];
115       const PetscReal y = points[p * dim + 1];
116       const PetscReal z = points[p * dim + 2];
117 
118       /* {2 y z, 0, 0} */
119       D[p * p_inc + 0 * b_inc + 0 * c_inc + 0] = 0.;
120       D[p * p_inc + 0 * b_inc + 0 * c_inc + 1] = 2. * z;
121       D[p * p_inc + 0 * b_inc + 0 * c_inc + 2] = 2. * y;
122       D[p * p_inc + 0 * b_inc + 1 * c_inc + 0] = 0.;
123       D[p * p_inc + 0 * b_inc + 1 * c_inc + 1] = 0.;
124       D[p * p_inc + 0 * b_inc + 1 * c_inc + 2] = 0.;
125       D[p * p_inc + 0 * b_inc + 2 * c_inc + 0] = 0.;
126       D[p * p_inc + 0 * b_inc + 2 * c_inc + 1] = 0.;
127       D[p * p_inc + 0 * b_inc + 2 * c_inc + 2] = 0.;
128       /* {0, 2 x z, 0} */
129       D[p * p_inc + 1 * b_inc + 0 * c_inc + 0] = 0.;
130       D[p * p_inc + 1 * b_inc + 0 * c_inc + 1] = 0.;
131       D[p * p_inc + 1 * b_inc + 0 * c_inc + 2] = 0.;
132       D[p * p_inc + 1 * b_inc + 1 * c_inc + 0] = 2. * z;
133       D[p * p_inc + 1 * b_inc + 1 * c_inc + 1] = 0.;
134       D[p * p_inc + 1 * b_inc + 1 * c_inc + 2] = 2. * x;
135       D[p * p_inc + 1 * b_inc + 2 * c_inc + 0] = 0.;
136       D[p * p_inc + 1 * b_inc + 2 * c_inc + 1] = 0.;
137       D[p * p_inc + 1 * b_inc + 2 * c_inc + 2] = 0.;
138       /* {0, 2 y z, -z^2} */
139       D[p * p_inc + 2 * b_inc + 0 * c_inc + 0] = 0.;
140       D[p * p_inc + 2 * b_inc + 0 * c_inc + 1] = 0.;
141       D[p * p_inc + 2 * b_inc + 0 * c_inc + 2] = 0.;
142       D[p * p_inc + 2 * b_inc + 1 * c_inc + 0] = 0.;
143       D[p * p_inc + 2 * b_inc + 1 * c_inc + 1] = 2. * z;
144       D[p * p_inc + 2 * b_inc + 1 * c_inc + 2] = 2. * y;
145       D[p * p_inc + 2 * b_inc + 2 * c_inc + 0] = 0.;
146       D[p * p_inc + 2 * b_inc + 2 * c_inc + 1] = 0.;
147       D[p * p_inc + 2 * b_inc + 2 * c_inc + 2] = -2. * z;
148       /* {2 x z, 0, -z^2} */
149       D[p * p_inc + 3 * b_inc + 0 * c_inc + 0] = 2. * z;
150       D[p * p_inc + 3 * b_inc + 0 * c_inc + 1] = 0.;
151       D[p * p_inc + 3 * b_inc + 0 * c_inc + 2] = 2. * x;
152       D[p * p_inc + 3 * b_inc + 1 * c_inc + 0] = 0.;
153       D[p * p_inc + 3 * b_inc + 1 * c_inc + 1] = 0.;
154       D[p * p_inc + 3 * b_inc + 1 * c_inc + 2] = 0.;
155       D[p * p_inc + 3 * b_inc + 2 * c_inc + 0] = 0.;
156       D[p * p_inc + 3 * b_inc + 2 * c_inc + 1] = 0.;
157       D[p * p_inc + 3 * b_inc + 2 * c_inc + 2] = -2. * z;
158       /* {x^2, x y, -3 x z} */
159       D[p * p_inc + 4 * b_inc + 0 * c_inc + 0] = 2. * x;
160       D[p * p_inc + 4 * b_inc + 0 * c_inc + 1] = 0.;
161       D[p * p_inc + 4 * b_inc + 0 * c_inc + 2] = 0.;
162       D[p * p_inc + 4 * b_inc + 1 * c_inc + 0] = y;
163       D[p * p_inc + 4 * b_inc + 1 * c_inc + 1] = x;
164       D[p * p_inc + 4 * b_inc + 1 * c_inc + 2] = 0.;
165       D[p * p_inc + 4 * b_inc + 2 * c_inc + 0] = -3. * z;
166       D[p * p_inc + 4 * b_inc + 2 * c_inc + 1] = 0.;
167       D[p * p_inc + 4 * b_inc + 2 * c_inc + 2] = -3. * x;
168       /* {x y, y^2, -3 y z} */
169       D[p * p_inc + 5 * b_inc + 0 * c_inc + 0] = y;
170       D[p * p_inc + 5 * b_inc + 0 * c_inc + 1] = x;
171       D[p * p_inc + 5 * b_inc + 0 * c_inc + 2] = 0.;
172       D[p * p_inc + 5 * b_inc + 1 * c_inc + 0] = 0.;
173       D[p * p_inc + 5 * b_inc + 1 * c_inc + 1] = 2. * y;
174       D[p * p_inc + 5 * b_inc + 1 * c_inc + 2] = 0.;
175       D[p * p_inc + 5 * b_inc + 2 * c_inc + 0] = 0.;
176       D[p * p_inc + 5 * b_inc + 2 * c_inc + 1] = -3. * z;
177       D[p * p_inc + 5 * b_inc + 2 * c_inc + 2] = -3. * y;
178     }
179   }
180   if (H) {
181     PetscInt p_inc = Nb * Nc * dim * dim;
182     PetscInt b_inc = Nc * dim * dim;
183     PetscInt c_inc = dim * dim;
184 
185     for (PetscInt p = 0; p < npoints; p++) {
186       /* {2 y z, 0, 0} */
187       D[p * p_inc + 0 * b_inc + 0 * c_inc + 0] = 0.;
188       D[p * p_inc + 0 * b_inc + 0 * c_inc + 1] = 0.;
189       D[p * p_inc + 0 * b_inc + 0 * c_inc + 2] = 0.;
190       D[p * p_inc + 0 * b_inc + 0 * c_inc + 3] = 0.;
191       D[p * p_inc + 0 * b_inc + 0 * c_inc + 4] = 0.;
192       D[p * p_inc + 0 * b_inc + 0 * c_inc + 5] = 2.;
193       D[p * p_inc + 0 * b_inc + 0 * c_inc + 6] = 0.;
194       D[p * p_inc + 0 * b_inc + 0 * c_inc + 7] = 2.;
195       D[p * p_inc + 0 * b_inc + 0 * c_inc + 8] = 0.;
196       D[p * p_inc + 0 * b_inc + 1 * c_inc + 0] = 0.;
197       D[p * p_inc + 0 * b_inc + 1 * c_inc + 1] = 0.;
198       D[p * p_inc + 0 * b_inc + 1 * c_inc + 2] = 0.;
199       D[p * p_inc + 0 * b_inc + 1 * c_inc + 3] = 0.;
200       D[p * p_inc + 0 * b_inc + 1 * c_inc + 4] = 0.;
201       D[p * p_inc + 0 * b_inc + 1 * c_inc + 5] = 0.;
202       D[p * p_inc + 0 * b_inc + 1 * c_inc + 6] = 0.;
203       D[p * p_inc + 0 * b_inc + 1 * c_inc + 7] = 0.;
204       D[p * p_inc + 0 * b_inc + 1 * c_inc + 8] = 0.;
205       D[p * p_inc + 0 * b_inc + 2 * c_inc + 0] = 0.;
206       D[p * p_inc + 0 * b_inc + 2 * c_inc + 1] = 0.;
207       D[p * p_inc + 0 * b_inc + 2 * c_inc + 2] = 0.;
208       D[p * p_inc + 0 * b_inc + 2 * c_inc + 3] = 0.;
209       D[p * p_inc + 0 * b_inc + 2 * c_inc + 4] = 0.;
210       D[p * p_inc + 0 * b_inc + 2 * c_inc + 5] = 0.;
211       D[p * p_inc + 0 * b_inc + 2 * c_inc + 6] = 0.;
212       D[p * p_inc + 0 * b_inc + 2 * c_inc + 7] = 0.;
213       D[p * p_inc + 0 * b_inc + 2 * c_inc + 8] = 0.;
214       /* {0, 2 x z, 0} */
215       D[p * p_inc + 1 * b_inc + 0 * c_inc + 0] = 0.;
216       D[p * p_inc + 1 * b_inc + 0 * c_inc + 1] = 0.;
217       D[p * p_inc + 1 * b_inc + 0 * c_inc + 2] = 0.;
218       D[p * p_inc + 1 * b_inc + 0 * c_inc + 3] = 0.;
219       D[p * p_inc + 1 * b_inc + 0 * c_inc + 4] = 0.;
220       D[p * p_inc + 1 * b_inc + 0 * c_inc + 5] = 0.;
221       D[p * p_inc + 1 * b_inc + 0 * c_inc + 6] = 0.;
222       D[p * p_inc + 1 * b_inc + 0 * c_inc + 7] = 0.;
223       D[p * p_inc + 1 * b_inc + 0 * c_inc + 8] = 0.;
224       D[p * p_inc + 1 * b_inc + 1 * c_inc + 0] = 0.;
225       D[p * p_inc + 1 * b_inc + 1 * c_inc + 1] = 0.;
226       D[p * p_inc + 1 * b_inc + 1 * c_inc + 2] = 2.;
227       D[p * p_inc + 1 * b_inc + 1 * c_inc + 3] = 0.;
228       D[p * p_inc + 1 * b_inc + 1 * c_inc + 4] = 0.;
229       D[p * p_inc + 1 * b_inc + 1 * c_inc + 5] = 0.;
230       D[p * p_inc + 1 * b_inc + 1 * c_inc + 6] = 2.;
231       D[p * p_inc + 1 * b_inc + 1 * c_inc + 7] = 0.;
232       D[p * p_inc + 1 * b_inc + 1 * c_inc + 8] = 0.;
233       D[p * p_inc + 1 * b_inc + 2 * c_inc + 0] = 0.;
234       D[p * p_inc + 1 * b_inc + 2 * c_inc + 1] = 0.;
235       D[p * p_inc + 1 * b_inc + 2 * c_inc + 2] = 0.;
236       D[p * p_inc + 1 * b_inc + 2 * c_inc + 3] = 0.;
237       D[p * p_inc + 1 * b_inc + 2 * c_inc + 4] = 0.;
238       D[p * p_inc + 1 * b_inc + 2 * c_inc + 5] = 0.;
239       D[p * p_inc + 1 * b_inc + 2 * c_inc + 6] = 0.;
240       D[p * p_inc + 1 * b_inc + 2 * c_inc + 7] = 0.;
241       D[p * p_inc + 1 * b_inc + 2 * c_inc + 8] = 0.;
242       /* {0, 2 y z, -z^2} */
243       D[p * p_inc + 2 * b_inc + 0 * c_inc + 0] = 0.;
244       D[p * p_inc + 2 * b_inc + 0 * c_inc + 1] = 0.;
245       D[p * p_inc + 2 * b_inc + 0 * c_inc + 2] = 0.;
246       D[p * p_inc + 2 * b_inc + 0 * c_inc + 3] = 0.;
247       D[p * p_inc + 2 * b_inc + 0 * c_inc + 4] = 0.;
248       D[p * p_inc + 2 * b_inc + 0 * c_inc + 5] = 0.;
249       D[p * p_inc + 2 * b_inc + 0 * c_inc + 6] = 0.;
250       D[p * p_inc + 2 * b_inc + 0 * c_inc + 7] = 0.;
251       D[p * p_inc + 2 * b_inc + 0 * c_inc + 8] = 0.;
252       D[p * p_inc + 2 * b_inc + 1 * c_inc + 0] = 0.;
253       D[p * p_inc + 2 * b_inc + 1 * c_inc + 1] = 0.;
254       D[p * p_inc + 2 * b_inc + 1 * c_inc + 2] = 0.;
255       D[p * p_inc + 2 * b_inc + 1 * c_inc + 3] = 0.;
256       D[p * p_inc + 2 * b_inc + 1 * c_inc + 4] = 0.;
257       D[p * p_inc + 2 * b_inc + 1 * c_inc + 5] = 2.;
258       D[p * p_inc + 2 * b_inc + 1 * c_inc + 6] = 0.;
259       D[p * p_inc + 2 * b_inc + 1 * c_inc + 7] = 2.;
260       D[p * p_inc + 2 * b_inc + 1 * c_inc + 8] = 0.;
261       D[p * p_inc + 2 * b_inc + 2 * c_inc + 0] = 0.;
262       D[p * p_inc + 2 * b_inc + 2 * c_inc + 1] = 0.;
263       D[p * p_inc + 2 * b_inc + 2 * c_inc + 2] = 0.;
264       D[p * p_inc + 2 * b_inc + 2 * c_inc + 3] = 0.;
265       D[p * p_inc + 2 * b_inc + 2 * c_inc + 4] = 0.;
266       D[p * p_inc + 2 * b_inc + 2 * c_inc + 5] = 0.;
267       D[p * p_inc + 2 * b_inc + 2 * c_inc + 6] = 0.;
268       D[p * p_inc + 2 * b_inc + 2 * c_inc + 7] = 0.;
269       D[p * p_inc + 2 * b_inc + 2 * c_inc + 8] = -2.;
270       /* {2 x z, 0, -z^2} */
271       D[p * p_inc + 3 * b_inc + 0 * c_inc + 0] = 0.;
272       D[p * p_inc + 3 * b_inc + 0 * c_inc + 1] = 0.;
273       D[p * p_inc + 3 * b_inc + 0 * c_inc + 2] = 2.;
274       D[p * p_inc + 3 * b_inc + 0 * c_inc + 3] = 0.;
275       D[p * p_inc + 3 * b_inc + 0 * c_inc + 4] = 0.;
276       D[p * p_inc + 3 * b_inc + 0 * c_inc + 5] = 0.;
277       D[p * p_inc + 3 * b_inc + 0 * c_inc + 6] = 2.;
278       D[p * p_inc + 3 * b_inc + 0 * c_inc + 7] = 0.;
279       D[p * p_inc + 3 * b_inc + 0 * c_inc + 8] = 0.;
280       D[p * p_inc + 3 * b_inc + 1 * c_inc + 0] = 0.;
281       D[p * p_inc + 3 * b_inc + 1 * c_inc + 1] = 0.;
282       D[p * p_inc + 3 * b_inc + 1 * c_inc + 2] = 0.;
283       D[p * p_inc + 3 * b_inc + 1 * c_inc + 3] = 0.;
284       D[p * p_inc + 3 * b_inc + 1 * c_inc + 4] = 0.;
285       D[p * p_inc + 3 * b_inc + 1 * c_inc + 5] = 0.;
286       D[p * p_inc + 3 * b_inc + 1 * c_inc + 6] = 0.;
287       D[p * p_inc + 3 * b_inc + 1 * c_inc + 7] = 0.;
288       D[p * p_inc + 3 * b_inc + 1 * c_inc + 8] = 0.;
289       D[p * p_inc + 3 * b_inc + 2 * c_inc + 0] = 0.;
290       D[p * p_inc + 3 * b_inc + 2 * c_inc + 1] = 0.;
291       D[p * p_inc + 3 * b_inc + 2 * c_inc + 2] = 0.;
292       D[p * p_inc + 3 * b_inc + 2 * c_inc + 3] = 0.;
293       D[p * p_inc + 3 * b_inc + 2 * c_inc + 4] = 0.;
294       D[p * p_inc + 3 * b_inc + 2 * c_inc + 5] = 0.;
295       D[p * p_inc + 3 * b_inc + 2 * c_inc + 6] = 0.;
296       D[p * p_inc + 3 * b_inc + 2 * c_inc + 7] = 0.;
297       D[p * p_inc + 3 * b_inc + 2 * c_inc + 8] = -2.;
298       /* {x^2, x y, -3 x z} */
299       D[p * p_inc + 4 * b_inc + 0 * c_inc + 0] = 2.;
300       D[p * p_inc + 4 * b_inc + 0 * c_inc + 1] = 0.;
301       D[p * p_inc + 4 * b_inc + 0 * c_inc + 2] = 0.;
302       D[p * p_inc + 4 * b_inc + 0 * c_inc + 3] = 0.;
303       D[p * p_inc + 4 * b_inc + 0 * c_inc + 4] = 0.;
304       D[p * p_inc + 4 * b_inc + 0 * c_inc + 5] = 0.;
305       D[p * p_inc + 4 * b_inc + 0 * c_inc + 6] = 0.;
306       D[p * p_inc + 4 * b_inc + 0 * c_inc + 7] = 0.;
307       D[p * p_inc + 4 * b_inc + 0 * c_inc + 8] = 0.;
308       D[p * p_inc + 4 * b_inc + 1 * c_inc + 0] = 0.;
309       D[p * p_inc + 4 * b_inc + 1 * c_inc + 1] = 1.;
310       D[p * p_inc + 4 * b_inc + 1 * c_inc + 2] = 0.;
311       D[p * p_inc + 4 * b_inc + 1 * c_inc + 3] = 1.;
312       D[p * p_inc + 4 * b_inc + 1 * c_inc + 4] = 0.;
313       D[p * p_inc + 4 * b_inc + 1 * c_inc + 5] = 0.;
314       D[p * p_inc + 4 * b_inc + 1 * c_inc + 6] = 0.;
315       D[p * p_inc + 4 * b_inc + 1 * c_inc + 7] = 0.;
316       D[p * p_inc + 4 * b_inc + 1 * c_inc + 8] = 0.;
317       D[p * p_inc + 4 * b_inc + 2 * c_inc + 0] = 0.;
318       D[p * p_inc + 4 * b_inc + 2 * c_inc + 1] = 0.;
319       D[p * p_inc + 4 * b_inc + 2 * c_inc + 2] = -3.;
320       D[p * p_inc + 4 * b_inc + 2 * c_inc + 3] = 0.;
321       D[p * p_inc + 4 * b_inc + 2 * c_inc + 4] = 0.;
322       D[p * p_inc + 4 * b_inc + 2 * c_inc + 5] = 0.;
323       D[p * p_inc + 4 * b_inc + 2 * c_inc + 6] = -3.;
324       D[p * p_inc + 4 * b_inc + 2 * c_inc + 7] = 0.;
325       D[p * p_inc + 4 * b_inc + 2 * c_inc + 8] = 0.;
326       /* {x y, y^2, -3 y z} */
327       D[p * p_inc + 4 * b_inc + 0 * c_inc + 0] = 2.;
328       D[p * p_inc + 4 * b_inc + 0 * c_inc + 1] = 0.;
329       D[p * p_inc + 4 * b_inc + 0 * c_inc + 2] = 0.;
330       D[p * p_inc + 4 * b_inc + 0 * c_inc + 3] = 0.;
331       D[p * p_inc + 4 * b_inc + 0 * c_inc + 4] = 0.;
332       D[p * p_inc + 4 * b_inc + 0 * c_inc + 5] = 0.;
333       D[p * p_inc + 4 * b_inc + 0 * c_inc + 6] = 0.;
334       D[p * p_inc + 4 * b_inc + 0 * c_inc + 7] = 0.;
335       D[p * p_inc + 4 * b_inc + 0 * c_inc + 8] = 0.;
336       D[p * p_inc + 4 * b_inc + 1 * c_inc + 0] = 0.;
337       D[p * p_inc + 4 * b_inc + 1 * c_inc + 1] = 1.;
338       D[p * p_inc + 4 * b_inc + 1 * c_inc + 2] = 0.;
339       D[p * p_inc + 4 * b_inc + 1 * c_inc + 3] = 1.;
340       D[p * p_inc + 4 * b_inc + 1 * c_inc + 4] = 0.;
341       D[p * p_inc + 4 * b_inc + 1 * c_inc + 5] = 0.;
342       D[p * p_inc + 4 * b_inc + 1 * c_inc + 6] = 0.;
343       D[p * p_inc + 4 * b_inc + 1 * c_inc + 7] = 0.;
344       D[p * p_inc + 4 * b_inc + 1 * c_inc + 8] = 0.;
345       D[p * p_inc + 4 * b_inc + 2 * c_inc + 0] = 0.;
346       D[p * p_inc + 4 * b_inc + 2 * c_inc + 1] = 0.;
347       D[p * p_inc + 4 * b_inc + 2 * c_inc + 2] = -3.;
348       D[p * p_inc + 4 * b_inc + 2 * c_inc + 3] = 0.;
349       D[p * p_inc + 4 * b_inc + 2 * c_inc + 4] = 0.;
350       D[p * p_inc + 4 * b_inc + 2 * c_inc + 5] = 0.;
351       D[p * p_inc + 4 * b_inc + 2 * c_inc + 6] = -3.;
352       D[p * p_inc + 4 * b_inc + 2 * c_inc + 7] = 0.;
353       D[p * p_inc + 4 * b_inc + 2 * c_inc + 8] = 0.;
354     }
355   }
356   PetscFunctionReturn(PETSC_SUCCESS);
357 }
358 
PetscSpaceGetHeightSubspace_WXY(PetscSpace sp,PetscInt height,PetscSpace * subsp)359 static PetscErrorCode PetscSpaceGetHeightSubspace_WXY(PetscSpace sp, PetscInt height, PetscSpace *subsp)
360 {
361   SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_SUP, "Do not know how to do this");
362 }
363 
PetscSpaceInitialize_WXY(PetscSpace sp)364 static PetscErrorCode PetscSpaceInitialize_WXY(PetscSpace sp)
365 {
366   PetscFunctionBegin;
367   sp->ops->setfromoptions    = PetscSpaceSetFromOptions_WXY;
368   sp->ops->setup             = PetscSpaceSetUp_WXY;
369   sp->ops->view              = PetscSpaceView_WXY;
370   sp->ops->destroy           = PetscSpaceDestroy_WXY;
371   sp->ops->getdimension      = PetscSpaceGetDimension_WXY;
372   sp->ops->evaluate          = PetscSpaceEvaluate_WXY;
373   sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_WXY;
374   PetscFunctionReturn(PETSC_SUCCESS);
375 }
376 
377 /*MC
378   PETSCSPACEWXY = "wxy" - A `PetscSpace` object that encapsulates the Wheeler-Xu-Yotov enrichments.
379 
380   Level: intermediate
381 
382   Note:
383 .vb
384   curl {{0, 0, y^2 z}, {x z^2, 0, 0}, {y z^2, 0, 0}, {0, -x z^2, 0}, {0, -3/2 x^2 z, -1/2 x^2 y}, {3/2 y^2 z, 0, 1/2 y^2 x}}
385   = {{2 y z, 0, 0}, {0, 2 x z, 0}, {0, 2 y z, -z^2}, {2 x z, 0, -z^2}, {x^2, x y, -3 x z}, {x y, y^2, -3 y z}}
386 .ve
387 
388 .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
389 M*/
390 
PetscSpaceCreate_WXY(PetscSpace sp)391 PETSC_EXTERN PetscErrorCode PetscSpaceCreate_WXY(PetscSpace sp)
392 {
393   PetscSpace_WXY *wxy;
394 
395   PetscFunctionBegin;
396   PetscValidHeaderSpecific(sp, PETSCSPACE_CLASSID, 1);
397   PetscCall(PetscNew(&wxy));
398   sp->data   = wxy;
399   sp->degree = 2;
400 
401   PetscCall(PetscSpaceInitialize_WXY(sp));
402   PetscFunctionReturn(PETSC_SUCCESS);
403 }
404