Actual source code: spacewxy.c
1: #include <petsc/private/petscfeimpl.h>
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: }
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: }
18: static PetscErrorCode PetscSpaceView_WXY(PetscSpace sp, PetscViewer viewer)
19: {
20: PetscBool iascii;
22: PetscFunctionBegin;
25: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
26: if (iascii) PetscCall(PetscSpacePolynomialView_Ascii(sp, viewer));
27: PetscFunctionReturn(PETSC_SUCCESS);
28: }
30: static PetscErrorCode PetscSpaceDestroy_WXY(PetscSpace sp)
31: {
32: PetscSpace_WXY *wxy = (PetscSpace_WXY *)sp->data;
34: PetscFunctionBegin;
35: PetscCall(PetscFree(wxy));
36: PetscFunctionReturn(PETSC_SUCCESS);
37: }
39: static PetscErrorCode PetscSpaceSetUp_WXY(PetscSpace sp)
40: {
41: PetscSpace_WXY *wxy = (PetscSpace_WXY *)sp->data;
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\n", sp->degree);
46: sp->maxDegree = sp->degree;
47: wxy->setupCalled = PETSC_TRUE;
48: PetscFunctionReturn(PETSC_SUCCESS);
49: }
51: static PetscErrorCode PetscSpaceGetDimension_WXY(PetscSpace sp, PetscInt *dim)
52: {
53: PetscFunctionBegin;
54: *dim = 6;
55: PetscFunctionReturn(PETSC_SUCCESS);
56: }
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;
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;
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];
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;
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];
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;
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: }
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: }
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: }
377: /*MC
378: PETSCSPACEWXY = "wxy" - A `PetscSpace` object that encapsulates the Wheeler-Xu-Yotov enrichments.
380: Level: intermediate
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
388: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
389: M*/
391: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_WXY(PetscSpace sp)
392: {
393: PetscSpace_WXY *wxy;
395: PetscFunctionBegin;
397: PetscCall(PetscNew(&wxy));
398: sp->data = wxy;
399: sp->degree = 2;
401: PetscCall(PetscSpaceInitialize_WXY(sp));
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }