Actual source code: spacewxy.c

  1: #include <petsc/private/petscfeimpl.h>

  3: static PetscErrorCode PetscSpaceSetFromOptions_WXY(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
  4: {
  5:   PetscOptionsHead(PetscOptionsObject,"PetscSpace WXY options");
  6:   PetscOptionsTail();
  7:   return 0;
  8: }

 10: static PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer v)
 11: {
 12:   PetscViewerASCIIPrintf(v, "WXY space of degree %" PetscInt_FMT "\n", sp->degree);
 13:   return 0;
 14: }

 16: static PetscErrorCode PetscSpaceView_WXY(PetscSpace sp, PetscViewer viewer)
 17: {
 18:   PetscBool      iascii;

 22:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
 23:   if (iascii) PetscSpacePolynomialView_Ascii(sp, viewer);
 24:   return 0;
 25: }

 27: static PetscErrorCode PetscSpaceDestroy_WXY(PetscSpace sp)
 28: {
 29:   PetscSpace_WXY *wxy = (PetscSpace_WXY *) sp->data;

 31:   PetscFree(wxy);
 32:   return 0;
 33: }

 35: static PetscErrorCode PetscSpaceSetUp_WXY(PetscSpace sp)
 36: {
 37:   PetscSpace_WXY *wxy = (PetscSpace_WXY *) sp->data;

 39:   if (wxy->setupCalled) return 0;
 41:   sp->maxDegree = sp->degree;
 42:   wxy->setupCalled = PETSC_TRUE;
 43:   return 0;
 44: }

 46: static PetscErrorCode PetscSpaceGetDimension_WXY(PetscSpace sp, PetscInt *dim)
 47: {
 48:   *dim = 6;
 49:   return 0;
 50: }

 52: static PetscErrorCode PetscSpaceEvaluate_WXY(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
 53: {
 54:   PetscSpace_WXY *wxy = (PetscSpace_WXY *) sp->data;
 55:   PetscInt        dim = sp->Nv;
 56:   PetscInt        Nb  = 6;
 57:   PetscInt        Nc  = 3;

 59:   if (!wxy->setupCalled) {
 60:     PetscSpaceSetUp(sp);
 61:     PetscSpaceEvaluate(sp, npoints, points, B, D, H);
 62:     return 0;
 63:   }
 65:   if (B) {
 66:     PetscInt p_inc = Nb*Nc;
 67:     PetscInt b_inc = Nc;
 68:     PetscInt c_inc = 1;

 70:     for (PetscInt p = 0; p < npoints; p++) {
 71:       const PetscReal x = points[p*dim+0];
 72:       const PetscReal y = points[p*dim+1];
 73:       const PetscReal z = points[p*dim+2];

 75:       /* {2 y z, 0, 0} */
 76:       B[p*p_inc + 0*b_inc + 0*c_inc] = 2.*y*z;
 77:       B[p*p_inc + 0*b_inc + 1*c_inc] = 0.;
 78:       B[p*p_inc + 0*b_inc + 2*c_inc] = 0.;
 79:       /* {0, 2 x z, 0} */
 80:       B[p*p_inc + 1*b_inc + 0*c_inc] = 0.;
 81:       B[p*p_inc + 1*b_inc + 1*c_inc] = 2.*x*z;
 82:       B[p*p_inc + 1*b_inc + 2*c_inc] = 0.;
 83:       /* {0, 2 y z, -z^2} */
 84:       B[p*p_inc + 2*b_inc + 0*c_inc] = 0.;
 85:       B[p*p_inc + 2*b_inc + 1*c_inc] = 2.*y*z;
 86:       B[p*p_inc + 2*b_inc + 2*c_inc] = -z*z;
 87:       /* {2 x z, 0, -z^2} */
 88:       B[p*p_inc + 3*b_inc + 0*c_inc] = 2.*x*z;
 89:       B[p*p_inc + 3*b_inc + 1*c_inc] = 0.;
 90:       B[p*p_inc + 3*b_inc + 2*c_inc] = -z*z;
 91:       /* {x^2, x y, -3 x z} */
 92:       B[p*p_inc + 4*b_inc + 0*c_inc] = x*x;
 93:       B[p*p_inc + 4*b_inc + 1*c_inc] = x*y;
 94:       B[p*p_inc + 4*b_inc + 2*c_inc] = -3.*x*z;
 95:       /* {x y, y^2, -3 y z} */
 96:       B[p*p_inc + 5*b_inc + 0*c_inc] = x*y;
 97:       B[p*p_inc + 5*b_inc + 1*c_inc] = y*y;
 98:       B[p*p_inc + 5*b_inc + 2*c_inc] = -3.*y*z;
 99:     }
100:   }
101:   if (D) {
102:     PetscInt p_inc = Nb*Nc*dim;
103:     PetscInt b_inc = Nc*dim;
104:     PetscInt c_inc = dim;

106:     for (PetscInt p = 0; p < npoints; p++) {
107:       const PetscReal x = points[p*dim+0];
108:       const PetscReal y = points[p*dim+1];
109:       const PetscReal z = points[p*dim+2];

111:       /* {2 y z, 0, 0} */
112:       D[p*p_inc + 0*b_inc + 0*c_inc + 0] = 0.;
113:       D[p*p_inc + 0*b_inc + 0*c_inc + 1] = 2.*z;
114:       D[p*p_inc + 0*b_inc + 0*c_inc + 2] = 2.*y;
115:       D[p*p_inc + 0*b_inc + 1*c_inc + 0] = 0.;
116:       D[p*p_inc + 0*b_inc + 1*c_inc + 1] = 0.;
117:       D[p*p_inc + 0*b_inc + 1*c_inc + 2] = 0.;
118:       D[p*p_inc + 0*b_inc + 2*c_inc + 0] = 0.;
119:       D[p*p_inc + 0*b_inc + 2*c_inc + 1] = 0.;
120:       D[p*p_inc + 0*b_inc + 2*c_inc + 2] = 0.;
121:       /* {0, 2 x z, 0} */
122:       D[p*p_inc + 1*b_inc + 0*c_inc + 0] = 0.;
123:       D[p*p_inc + 1*b_inc + 0*c_inc + 1] = 0.;
124:       D[p*p_inc + 1*b_inc + 0*c_inc + 2] = 0.;
125:       D[p*p_inc + 1*b_inc + 1*c_inc + 0] = 2.*z;
126:       D[p*p_inc + 1*b_inc + 1*c_inc + 1] = 0.;
127:       D[p*p_inc + 1*b_inc + 1*c_inc + 2] = 2.*x;
128:       D[p*p_inc + 1*b_inc + 2*c_inc + 0] = 0.;
129:       D[p*p_inc + 1*b_inc + 2*c_inc + 1] = 0.;
130:       D[p*p_inc + 1*b_inc + 2*c_inc + 2] = 0.;
131:       /* {0, 2 y z, -z^2} */
132:       D[p*p_inc + 2*b_inc + 0*c_inc + 0] = 0.;
133:       D[p*p_inc + 2*b_inc + 0*c_inc + 1] = 0.;
134:       D[p*p_inc + 2*b_inc + 0*c_inc + 2] = 0.;
135:       D[p*p_inc + 2*b_inc + 1*c_inc + 0] = 0.;
136:       D[p*p_inc + 2*b_inc + 1*c_inc + 1] = 2.*z;
137:       D[p*p_inc + 2*b_inc + 1*c_inc + 2] = 2.*y;
138:       D[p*p_inc + 2*b_inc + 2*c_inc + 0] = 0.;
139:       D[p*p_inc + 2*b_inc + 2*c_inc + 1] = 0.;
140:       D[p*p_inc + 2*b_inc + 2*c_inc + 2] = -2.*z;
141:       /* {2 x z, 0, -z^2} */
142:       D[p*p_inc + 3*b_inc + 0*c_inc + 0] = 2.*z;
143:       D[p*p_inc + 3*b_inc + 0*c_inc + 1] = 0.;
144:       D[p*p_inc + 3*b_inc + 0*c_inc + 2] = 2.*x;
145:       D[p*p_inc + 3*b_inc + 1*c_inc + 0] = 0.;
146:       D[p*p_inc + 3*b_inc + 1*c_inc + 1] = 0.;
147:       D[p*p_inc + 3*b_inc + 1*c_inc + 2] = 0.;
148:       D[p*p_inc + 3*b_inc + 2*c_inc + 0] = 0.;
149:       D[p*p_inc + 3*b_inc + 2*c_inc + 1] = 0.;
150:       D[p*p_inc + 3*b_inc + 2*c_inc + 2] = -2.*z;
151:       /* {x^2, x y, -3 x z} */
152:       D[p*p_inc + 4*b_inc + 0*c_inc + 0] = 2.*x;
153:       D[p*p_inc + 4*b_inc + 0*c_inc + 1] = 0.;
154:       D[p*p_inc + 4*b_inc + 0*c_inc + 2] = 0.;
155:       D[p*p_inc + 4*b_inc + 1*c_inc + 0] = y;
156:       D[p*p_inc + 4*b_inc + 1*c_inc + 1] = x;
157:       D[p*p_inc + 4*b_inc + 1*c_inc + 2] = 0.;
158:       D[p*p_inc + 4*b_inc + 2*c_inc + 0] = -3.*z;
159:       D[p*p_inc + 4*b_inc + 2*c_inc + 1] = 0.;
160:       D[p*p_inc + 4*b_inc + 2*c_inc + 2] = -3.*x;
161:       /* {x y, y^2, -3 y z} */
162:       D[p*p_inc + 5*b_inc + 0*c_inc + 0] = y;
163:       D[p*p_inc + 5*b_inc + 0*c_inc + 1] = x;
164:       D[p*p_inc + 5*b_inc + 0*c_inc + 2] = 0.;
165:       D[p*p_inc + 5*b_inc + 1*c_inc + 0] = 0.;
166:       D[p*p_inc + 5*b_inc + 1*c_inc + 1] = 2.*y;
167:       D[p*p_inc + 5*b_inc + 1*c_inc + 2] = 0.;
168:       D[p*p_inc + 5*b_inc + 2*c_inc + 0] = 0.;
169:       D[p*p_inc + 5*b_inc + 2*c_inc + 1] = -3.*z;
170:       D[p*p_inc + 5*b_inc + 2*c_inc + 2] = -3.*y;
171:     }
172:   }
173:   if (H) {
174:     PetscInt p_inc = Nb*Nc*dim*dim;
175:     PetscInt b_inc = Nc*dim*dim;
176:     PetscInt c_inc = dim*dim;

178:     for (PetscInt p = 0; p < npoints; p++) {
179:       /* {2 y z, 0, 0} */
180:       D[p*p_inc + 0*b_inc + 0*c_inc + 0] = 0.;
181:       D[p*p_inc + 0*b_inc + 0*c_inc + 1] = 0.;
182:       D[p*p_inc + 0*b_inc + 0*c_inc + 2] = 0.;
183:       D[p*p_inc + 0*b_inc + 0*c_inc + 3] = 0.;
184:       D[p*p_inc + 0*b_inc + 0*c_inc + 4] = 0.;
185:       D[p*p_inc + 0*b_inc + 0*c_inc + 5] = 2.;
186:       D[p*p_inc + 0*b_inc + 0*c_inc + 6] = 0.;
187:       D[p*p_inc + 0*b_inc + 0*c_inc + 7] = 2.;
188:       D[p*p_inc + 0*b_inc + 0*c_inc + 8] = 0.;
189:       D[p*p_inc + 0*b_inc + 1*c_inc + 0] = 0.;
190:       D[p*p_inc + 0*b_inc + 1*c_inc + 1] = 0.;
191:       D[p*p_inc + 0*b_inc + 1*c_inc + 2] = 0.;
192:       D[p*p_inc + 0*b_inc + 1*c_inc + 3] = 0.;
193:       D[p*p_inc + 0*b_inc + 1*c_inc + 4] = 0.;
194:       D[p*p_inc + 0*b_inc + 1*c_inc + 5] = 0.;
195:       D[p*p_inc + 0*b_inc + 1*c_inc + 6] = 0.;
196:       D[p*p_inc + 0*b_inc + 1*c_inc + 7] = 0.;
197:       D[p*p_inc + 0*b_inc + 1*c_inc + 8] = 0.;
198:       D[p*p_inc + 0*b_inc + 2*c_inc + 0] = 0.;
199:       D[p*p_inc + 0*b_inc + 2*c_inc + 1] = 0.;
200:       D[p*p_inc + 0*b_inc + 2*c_inc + 2] = 0.;
201:       D[p*p_inc + 0*b_inc + 2*c_inc + 3] = 0.;
202:       D[p*p_inc + 0*b_inc + 2*c_inc + 4] = 0.;
203:       D[p*p_inc + 0*b_inc + 2*c_inc + 5] = 0.;
204:       D[p*p_inc + 0*b_inc + 2*c_inc + 6] = 0.;
205:       D[p*p_inc + 0*b_inc + 2*c_inc + 7] = 0.;
206:       D[p*p_inc + 0*b_inc + 2*c_inc + 8] = 0.;
207:       /* {0, 2 x z, 0} */
208:       D[p*p_inc + 1*b_inc + 0*c_inc + 0] = 0.;
209:       D[p*p_inc + 1*b_inc + 0*c_inc + 1] = 0.;
210:       D[p*p_inc + 1*b_inc + 0*c_inc + 2] = 0.;
211:       D[p*p_inc + 1*b_inc + 0*c_inc + 3] = 0.;
212:       D[p*p_inc + 1*b_inc + 0*c_inc + 4] = 0.;
213:       D[p*p_inc + 1*b_inc + 0*c_inc + 5] = 0.;
214:       D[p*p_inc + 1*b_inc + 0*c_inc + 6] = 0.;
215:       D[p*p_inc + 1*b_inc + 0*c_inc + 7] = 0.;
216:       D[p*p_inc + 1*b_inc + 0*c_inc + 8] = 0.;
217:       D[p*p_inc + 1*b_inc + 1*c_inc + 0] = 0.;
218:       D[p*p_inc + 1*b_inc + 1*c_inc + 1] = 0.;
219:       D[p*p_inc + 1*b_inc + 1*c_inc + 2] = 2.;
220:       D[p*p_inc + 1*b_inc + 1*c_inc + 3] = 0.;
221:       D[p*p_inc + 1*b_inc + 1*c_inc + 4] = 0.;
222:       D[p*p_inc + 1*b_inc + 1*c_inc + 5] = 0.;
223:       D[p*p_inc + 1*b_inc + 1*c_inc + 6] = 2.;
224:       D[p*p_inc + 1*b_inc + 1*c_inc + 7] = 0.;
225:       D[p*p_inc + 1*b_inc + 1*c_inc + 8] = 0.;
226:       D[p*p_inc + 1*b_inc + 2*c_inc + 0] = 0.;
227:       D[p*p_inc + 1*b_inc + 2*c_inc + 1] = 0.;
228:       D[p*p_inc + 1*b_inc + 2*c_inc + 2] = 0.;
229:       D[p*p_inc + 1*b_inc + 2*c_inc + 3] = 0.;
230:       D[p*p_inc + 1*b_inc + 2*c_inc + 4] = 0.;
231:       D[p*p_inc + 1*b_inc + 2*c_inc + 5] = 0.;
232:       D[p*p_inc + 1*b_inc + 2*c_inc + 6] = 0.;
233:       D[p*p_inc + 1*b_inc + 2*c_inc + 7] = 0.;
234:       D[p*p_inc + 1*b_inc + 2*c_inc + 8] = 0.;
235:       /* {0, 2 y z, -z^2} */
236:       D[p*p_inc + 2*b_inc + 0*c_inc + 0] = 0.;
237:       D[p*p_inc + 2*b_inc + 0*c_inc + 1] = 0.;
238:       D[p*p_inc + 2*b_inc + 0*c_inc + 2] = 0.;
239:       D[p*p_inc + 2*b_inc + 0*c_inc + 3] = 0.;
240:       D[p*p_inc + 2*b_inc + 0*c_inc + 4] = 0.;
241:       D[p*p_inc + 2*b_inc + 0*c_inc + 5] = 0.;
242:       D[p*p_inc + 2*b_inc + 0*c_inc + 6] = 0.;
243:       D[p*p_inc + 2*b_inc + 0*c_inc + 7] = 0.;
244:       D[p*p_inc + 2*b_inc + 0*c_inc + 8] = 0.;
245:       D[p*p_inc + 2*b_inc + 1*c_inc + 0] = 0.;
246:       D[p*p_inc + 2*b_inc + 1*c_inc + 1] = 0.;
247:       D[p*p_inc + 2*b_inc + 1*c_inc + 2] = 0.;
248:       D[p*p_inc + 2*b_inc + 1*c_inc + 3] = 0.;
249:       D[p*p_inc + 2*b_inc + 1*c_inc + 4] = 0.;
250:       D[p*p_inc + 2*b_inc + 1*c_inc + 5] = 2.;
251:       D[p*p_inc + 2*b_inc + 1*c_inc + 6] = 0.;
252:       D[p*p_inc + 2*b_inc + 1*c_inc + 7] = 2.;
253:       D[p*p_inc + 2*b_inc + 1*c_inc + 8] = 0.;
254:       D[p*p_inc + 2*b_inc + 2*c_inc + 0] = 0.;
255:       D[p*p_inc + 2*b_inc + 2*c_inc + 1] = 0.;
256:       D[p*p_inc + 2*b_inc + 2*c_inc + 2] = 0.;
257:       D[p*p_inc + 2*b_inc + 2*c_inc + 3] = 0.;
258:       D[p*p_inc + 2*b_inc + 2*c_inc + 4] = 0.;
259:       D[p*p_inc + 2*b_inc + 2*c_inc + 5] = 0.;
260:       D[p*p_inc + 2*b_inc + 2*c_inc + 6] = 0.;
261:       D[p*p_inc + 2*b_inc + 2*c_inc + 7] = 0.;
262:       D[p*p_inc + 2*b_inc + 2*c_inc + 8] = -2.;
263:       /* {2 x z, 0, -z^2} */
264:       D[p*p_inc + 3*b_inc + 0*c_inc + 0] = 0.;
265:       D[p*p_inc + 3*b_inc + 0*c_inc + 1] = 0.;
266:       D[p*p_inc + 3*b_inc + 0*c_inc + 2] = 2.;
267:       D[p*p_inc + 3*b_inc + 0*c_inc + 3] = 0.;
268:       D[p*p_inc + 3*b_inc + 0*c_inc + 4] = 0.;
269:       D[p*p_inc + 3*b_inc + 0*c_inc + 5] = 0.;
270:       D[p*p_inc + 3*b_inc + 0*c_inc + 6] = 2.;
271:       D[p*p_inc + 3*b_inc + 0*c_inc + 7] = 0.;
272:       D[p*p_inc + 3*b_inc + 0*c_inc + 8] = 0.;
273:       D[p*p_inc + 3*b_inc + 1*c_inc + 0] = 0.;
274:       D[p*p_inc + 3*b_inc + 1*c_inc + 1] = 0.;
275:       D[p*p_inc + 3*b_inc + 1*c_inc + 2] = 0.;
276:       D[p*p_inc + 3*b_inc + 1*c_inc + 3] = 0.;
277:       D[p*p_inc + 3*b_inc + 1*c_inc + 4] = 0.;
278:       D[p*p_inc + 3*b_inc + 1*c_inc + 5] = 0.;
279:       D[p*p_inc + 3*b_inc + 1*c_inc + 6] = 0.;
280:       D[p*p_inc + 3*b_inc + 1*c_inc + 7] = 0.;
281:       D[p*p_inc + 3*b_inc + 1*c_inc + 8] = 0.;
282:       D[p*p_inc + 3*b_inc + 2*c_inc + 0] = 0.;
283:       D[p*p_inc + 3*b_inc + 2*c_inc + 1] = 0.;
284:       D[p*p_inc + 3*b_inc + 2*c_inc + 2] = 0.;
285:       D[p*p_inc + 3*b_inc + 2*c_inc + 3] = 0.;
286:       D[p*p_inc + 3*b_inc + 2*c_inc + 4] = 0.;
287:       D[p*p_inc + 3*b_inc + 2*c_inc + 5] = 0.;
288:       D[p*p_inc + 3*b_inc + 2*c_inc + 6] = 0.;
289:       D[p*p_inc + 3*b_inc + 2*c_inc + 7] = 0.;
290:       D[p*p_inc + 3*b_inc + 2*c_inc + 8] = -2.;
291:       /* {x^2, x y, -3 x z} */
292:       D[p*p_inc + 4*b_inc + 0*c_inc + 0] = 2.;
293:       D[p*p_inc + 4*b_inc + 0*c_inc + 1] = 0.;
294:       D[p*p_inc + 4*b_inc + 0*c_inc + 2] = 0.;
295:       D[p*p_inc + 4*b_inc + 0*c_inc + 3] = 0.;
296:       D[p*p_inc + 4*b_inc + 0*c_inc + 4] = 0.;
297:       D[p*p_inc + 4*b_inc + 0*c_inc + 5] = 0.;
298:       D[p*p_inc + 4*b_inc + 0*c_inc + 6] = 0.;
299:       D[p*p_inc + 4*b_inc + 0*c_inc + 7] = 0.;
300:       D[p*p_inc + 4*b_inc + 0*c_inc + 8] = 0.;
301:       D[p*p_inc + 4*b_inc + 1*c_inc + 0] = 0.;
302:       D[p*p_inc + 4*b_inc + 1*c_inc + 1] = 1.;
303:       D[p*p_inc + 4*b_inc + 1*c_inc + 2] = 0.;
304:       D[p*p_inc + 4*b_inc + 1*c_inc + 3] = 1.;
305:       D[p*p_inc + 4*b_inc + 1*c_inc + 4] = 0.;
306:       D[p*p_inc + 4*b_inc + 1*c_inc + 5] = 0.;
307:       D[p*p_inc + 4*b_inc + 1*c_inc + 6] = 0.;
308:       D[p*p_inc + 4*b_inc + 1*c_inc + 7] = 0.;
309:       D[p*p_inc + 4*b_inc + 1*c_inc + 8] = 0.;
310:       D[p*p_inc + 4*b_inc + 2*c_inc + 0] = 0.;
311:       D[p*p_inc + 4*b_inc + 2*c_inc + 1] = 0.;
312:       D[p*p_inc + 4*b_inc + 2*c_inc + 2] = -3.;
313:       D[p*p_inc + 4*b_inc + 2*c_inc + 3] = 0.;
314:       D[p*p_inc + 4*b_inc + 2*c_inc + 4] = 0.;
315:       D[p*p_inc + 4*b_inc + 2*c_inc + 5] = 0.;
316:       D[p*p_inc + 4*b_inc + 2*c_inc + 6] = -3.;
317:       D[p*p_inc + 4*b_inc + 2*c_inc + 7] = 0.;
318:       D[p*p_inc + 4*b_inc + 2*c_inc + 8] = 0.;
319:       /* {x y, y^2, -3 y z} */
320:       D[p*p_inc + 4*b_inc + 0*c_inc + 0] = 2.;
321:       D[p*p_inc + 4*b_inc + 0*c_inc + 1] = 0.;
322:       D[p*p_inc + 4*b_inc + 0*c_inc + 2] = 0.;
323:       D[p*p_inc + 4*b_inc + 0*c_inc + 3] = 0.;
324:       D[p*p_inc + 4*b_inc + 0*c_inc + 4] = 0.;
325:       D[p*p_inc + 4*b_inc + 0*c_inc + 5] = 0.;
326:       D[p*p_inc + 4*b_inc + 0*c_inc + 6] = 0.;
327:       D[p*p_inc + 4*b_inc + 0*c_inc + 7] = 0.;
328:       D[p*p_inc + 4*b_inc + 0*c_inc + 8] = 0.;
329:       D[p*p_inc + 4*b_inc + 1*c_inc + 0] = 0.;
330:       D[p*p_inc + 4*b_inc + 1*c_inc + 1] = 1.;
331:       D[p*p_inc + 4*b_inc + 1*c_inc + 2] = 0.;
332:       D[p*p_inc + 4*b_inc + 1*c_inc + 3] = 1.;
333:       D[p*p_inc + 4*b_inc + 1*c_inc + 4] = 0.;
334:       D[p*p_inc + 4*b_inc + 1*c_inc + 5] = 0.;
335:       D[p*p_inc + 4*b_inc + 1*c_inc + 6] = 0.;
336:       D[p*p_inc + 4*b_inc + 1*c_inc + 7] = 0.;
337:       D[p*p_inc + 4*b_inc + 1*c_inc + 8] = 0.;
338:       D[p*p_inc + 4*b_inc + 2*c_inc + 0] = 0.;
339:       D[p*p_inc + 4*b_inc + 2*c_inc + 1] = 0.;
340:       D[p*p_inc + 4*b_inc + 2*c_inc + 2] = -3.;
341:       D[p*p_inc + 4*b_inc + 2*c_inc + 3] = 0.;
342:       D[p*p_inc + 4*b_inc + 2*c_inc + 4] = 0.;
343:       D[p*p_inc + 4*b_inc + 2*c_inc + 5] = 0.;
344:       D[p*p_inc + 4*b_inc + 2*c_inc + 6] = -3.;
345:       D[p*p_inc + 4*b_inc + 2*c_inc + 7] = 0.;
346:       D[p*p_inc + 4*b_inc + 2*c_inc + 8] = 0.;
347:     }
348:   }
349:   return 0;
350: }

352: static PetscErrorCode PetscSpaceGetHeightSubspace_WXY(PetscSpace sp, PetscInt height, PetscSpace *subsp)
353: {
354:   SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Do not know how to do this");
355: }

357: static PetscErrorCode PetscSpaceInitialize_WXY(PetscSpace sp)
358: {
359:   sp->ops->setfromoptions    = PetscSpaceSetFromOptions_WXY;
360:   sp->ops->setup             = PetscSpaceSetUp_WXY;
361:   sp->ops->view              = PetscSpaceView_WXY;
362:   sp->ops->destroy           = PetscSpaceDestroy_WXY;
363:   sp->ops->getdimension      = PetscSpaceGetDimension_WXY;
364:   sp->ops->evaluate          = PetscSpaceEvaluate_WXY;
365:   sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_WXY;
366:   return 0;
367: }

369: /*MC
370:   PETSCSPACEWXY = "wxy" - A PetscSpace object that encapsulates the Wheeler-Xu-Yotov enrichments.
371: $ 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}}
372: $ = {{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}}

374:   Level: intermediate

376: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
377: M*/

379: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_WXY(PetscSpace sp)
380: {
381:   PetscSpace_WXY *wxy;

384:   PetscNewLog(sp,&wxy);
385:   sp->data = wxy;
386:   sp->degree = 2;

388:   PetscSpaceInitialize_WXY(sp);
389:   return 0;
390: }