Actual source code: spacetensor.c
1: #include <petsc/private/petscfeimpl.h>
3: static PetscErrorCode PetscSpaceTensorCreateSubspace(PetscSpace space, PetscInt Nvs, PetscInt Ncs, PetscSpace *subspace)
4: {
5: PetscInt degree;
6: const char *prefix;
7: const char *name;
8: char subname[PETSC_MAX_PATH_LEN];
10: PetscSpaceGetDegree(space, °ree, NULL);
11: PetscObjectGetOptionsPrefix((PetscObject)space, &prefix);
12: PetscSpaceCreate(PetscObjectComm((PetscObject)space), subspace);
13: PetscSpaceSetType(*subspace, PETSCSPACEPOLYNOMIAL);
14: PetscSpaceSetNumVariables(*subspace, Nvs);
15: PetscSpaceSetNumComponents(*subspace, Ncs);
16: PetscSpaceSetDegree(*subspace, degree, PETSC_DETERMINE);
17: PetscObjectSetOptionsPrefix((PetscObject)*subspace, prefix);
18: PetscObjectAppendOptionsPrefix((PetscObject)*subspace, "tensorcomp_");
19: if (((PetscObject) space)->name) {
20: PetscObjectGetName((PetscObject)space, &name);
21: PetscSNPrintf(subname, PETSC_MAX_PATH_LEN-1, "%s tensor component", name);
22: PetscObjectSetName((PetscObject)*subspace, subname);
23: } else {
24: PetscObjectSetName((PetscObject)*subspace, "tensor component");
25: }
26: return 0;
27: }
29: static PetscErrorCode PetscSpaceSetFromOptions_Tensor(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
30: {
31: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
32: PetscInt Ns, Nc, i, Nv, deg;
33: PetscBool uniform = PETSC_TRUE;
35: PetscSpaceGetNumVariables(sp, &Nv);
36: if (!Nv) return 0;
37: PetscSpaceGetNumComponents(sp, &Nc);
38: PetscSpaceTensorGetNumSubspaces(sp, &Ns);
39: PetscSpaceGetDegree(sp, °, NULL);
40: if (Ns > 1) {
41: PetscSpace s0;
43: PetscSpaceTensorGetSubspace(sp, 0, &s0);
44: for (i = 1; i < Ns; i++) {
45: PetscSpace si;
47: PetscSpaceTensorGetSubspace(sp, i, &si);
48: if (si != s0) {uniform = PETSC_FALSE; break;}
49: }
50: }
51: Ns = (Ns == PETSC_DEFAULT) ? PetscMax(Nv,1) : Ns;
52: PetscOptionsHead(PetscOptionsObject,"PetscSpace tensor options");
53: PetscOptionsBoundedInt("-petscspace_tensor_spaces", "The number of subspaces", "PetscSpaceTensorSetNumSubspaces", Ns, &Ns, NULL,0);
54: PetscOptionsBool("-petscspace_tensor_uniform", "Subspaces are identical", "PetscSpaceTensorSetFromOptions", uniform, &uniform, NULL);
55: PetscOptionsTail();
58: if (Ns != tens->numTensSpaces) PetscSpaceTensorSetNumSubspaces(sp, Ns);
59: if (uniform) {
60: PetscInt Nvs = Nv / Ns;
61: PetscInt Ncs;
62: PetscSpace subspace;
65: Ncs = (PetscInt) PetscPowReal((PetscReal) Nc, 1./Ns);
67: PetscSpaceTensorGetSubspace(sp, 0, &subspace);
68: if (!subspace) PetscSpaceTensorCreateSubspace(sp, Nvs, Ncs, &subspace);
69: else PetscObjectReference((PetscObject)subspace);
70: PetscSpaceSetFromOptions(subspace);
71: for (i = 0; i < Ns; i++) PetscSpaceTensorSetSubspace(sp, i, subspace);
72: PetscSpaceDestroy(&subspace);
73: } else {
74: for (i = 0; i < Ns; i++) {
75: PetscSpace subspace;
77: PetscSpaceTensorGetSubspace(sp, i, &subspace);
78: if (!subspace) {
79: char tprefix[128];
81: PetscSpaceTensorCreateSubspace(sp, 1, 1, &subspace);
82: PetscSNPrintf(tprefix, 128, "%d_",(int)i);
83: PetscObjectAppendOptionsPrefix((PetscObject)subspace, tprefix);
84: } else {
85: PetscObjectReference((PetscObject)subspace);
86: }
87: PetscSpaceSetFromOptions(subspace);
88: PetscSpaceTensorSetSubspace(sp, i, subspace);
89: PetscSpaceDestroy(&subspace);
90: }
91: }
92: return 0;
93: }
95: static PetscErrorCode PetscSpaceTensorView_Ascii(PetscSpace sp, PetscViewer v)
96: {
97: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
98: PetscBool uniform = PETSC_TRUE;
99: PetscInt Ns = tens->numTensSpaces, i, n;
101: for (i = 1; i < Ns; i++) {
102: if (tens->tensspaces[i] != tens->tensspaces[0]) {uniform = PETSC_FALSE; break;}
103: }
104: if (uniform) PetscViewerASCIIPrintf(v, "Tensor space of %D subspaces (all identical)\n", Ns);
105: else PetscViewerASCIIPrintf(v, "Tensor space of %D subspaces\n", Ns);
106: n = uniform ? 1 : Ns;
107: for (i = 0; i < n; i++) {
108: PetscViewerASCIIPushTab(v);
109: PetscSpaceView(tens->tensspaces[i], v);
110: PetscViewerASCIIPopTab(v);
111: }
112: return 0;
113: }
115: static PetscErrorCode PetscSpaceView_Tensor(PetscSpace sp, PetscViewer viewer)
116: {
117: PetscBool iascii;
119: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
120: if (iascii) PetscSpaceTensorView_Ascii(sp, viewer);
121: return 0;
122: }
124: static PetscErrorCode PetscSpaceSetUp_Tensor(PetscSpace sp)
125: {
126: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
127: PetscInt Nc, Nv, Ns;
128: PetscBool uniform = PETSC_TRUE;
129: PetscInt deg, maxDeg;
130: PetscInt Ncprod;
132: if (tens->setupCalled) return 0;
133: PetscSpaceGetNumVariables(sp, &Nv);
134: PetscSpaceGetNumComponents(sp, &Nc);
135: PetscSpaceTensorGetNumSubspaces(sp, &Ns);
136: if (Ns == PETSC_DEFAULT) {
137: Ns = Nv;
138: PetscSpaceTensorSetNumSubspaces(sp, Ns);
139: }
140: if (!Ns) {
141: SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zero subspaces");
142: } else {
143: PetscSpace s0;
146: PetscSpaceTensorGetSubspace(sp, 0, &s0);
147: for (PetscInt i = 1; i < Ns; i++) {
148: PetscSpace si;
150: PetscSpaceTensorGetSubspace(sp, i, &si);
151: if (si != s0) {uniform = PETSC_FALSE; break;}
152: }
153: if (uniform) {
154: PetscInt Nvs = Nv / Ns;
155: PetscInt Ncs;
158: Ncs = (PetscInt) (PetscPowReal((PetscReal) Nc, 1./Ns));
160: if (!s0) PetscSpaceTensorCreateSubspace(sp, Nvs, Ncs, &s0);
161: else PetscObjectReference((PetscObject) s0);
162: PetscSpaceSetUp(s0);
163: for (PetscInt i = 0; i < Ns; i++) PetscSpaceTensorSetSubspace(sp, i, s0);
164: PetscSpaceDestroy(&s0);
165: Ncprod = PetscPowInt(Ncs, Ns);
166: } else {
167: PetscInt Nvsum = 0;
169: Ncprod = 1;
170: for (PetscInt i = 0 ; i < Ns; i++) {
171: PetscInt Nvs, Ncs;
172: PetscSpace si;
174: PetscSpaceTensorGetSubspace(sp, i, &si);
175: if (!si) PetscSpaceTensorCreateSubspace(sp, 1, 1, &si);
176: else PetscObjectReference((PetscObject) si);
177: PetscSpaceSetUp(si);
178: PetscSpaceTensorSetSubspace(sp, i, si);
179: PetscSpaceGetNumVariables(si, &Nvs);
180: PetscSpaceGetNumComponents(si, &Ncs);
181: PetscSpaceDestroy(&si);
183: Nvsum += Nvs;
184: Ncprod *= Ncs;
185: }
189: }
190: if (Ncprod != Nc) {
191: PetscInt Ncopies = Nc / Ncprod;
192: PetscInt Nv = sp->Nv;
193: const char *prefix;
194: const char *name;
195: char subname[PETSC_MAX_PATH_LEN];
196: PetscSpace subsp;
198: PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp);
199: PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix);
200: PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix);
201: PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_");
202: if (((PetscObject)sp)->name) {
203: PetscObjectGetName((PetscObject)sp, &name);
204: PetscSNPrintf(subname, PETSC_MAX_PATH_LEN-1, "%s sum component", name);
205: PetscObjectSetName((PetscObject)subsp, subname);
206: } else {
207: PetscObjectSetName((PetscObject)subsp, "sum component");
208: }
209: PetscSpaceSetType(subsp, PETSCSPACETENSOR);
210: PetscSpaceSetNumVariables(subsp, Nv);
211: PetscSpaceSetNumComponents(subsp, Ncprod);
212: PetscSpaceTensorSetNumSubspaces(subsp, Ns);
213: for (PetscInt i = 0; i < Ns; i++) {
214: PetscSpace ssp;
216: PetscSpaceTensorGetSubspace(sp, i, &ssp);
217: PetscSpaceTensorSetSubspace(subsp, i, ssp);
218: }
219: PetscSpaceSetUp(subsp);
220: PetscSpaceSetType(sp, PETSCSPACESUM);
221: PetscSpaceSumSetNumSubspaces(sp, Ncopies);
222: for (PetscInt i = 0; i < Ncopies; i++) {
223: PetscSpaceSumSetSubspace(sp, i, subsp);
224: }
225: PetscSpaceDestroy(&subsp);
226: PetscSpaceSetUp(sp);
227: return 0;
228: }
229: }
230: deg = PETSC_MAX_INT;
231: maxDeg = 0;
232: for (PetscInt i = 0; i < Ns; i++) {
233: PetscSpace si;
234: PetscInt iDeg, iMaxDeg;
236: PetscSpaceTensorGetSubspace(sp, i, &si);
237: PetscSpaceGetDegree(si, &iDeg, &iMaxDeg);
238: deg = PetscMin(deg, iDeg);
239: maxDeg += iMaxDeg;
240: }
241: sp->degree = deg;
242: sp->maxDegree = maxDeg;
243: tens->uniform = uniform;
244: tens->setupCalled = PETSC_TRUE;
245: return 0;
246: }
248: static PetscErrorCode PetscSpaceDestroy_Tensor(PetscSpace sp)
249: {
250: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
251: PetscInt Ns, i;
253: Ns = tens->numTensSpaces;
254: if (tens->heightsubspaces) {
255: PetscInt d;
257: /* sp->Nv is the spatial dimension, so it is equal to the number
258: * of subspaces on higher co-dimension points */
259: for (d = 0; d < sp->Nv; ++d) {
260: PetscSpaceDestroy(&tens->heightsubspaces[d]);
261: }
262: }
263: PetscFree(tens->heightsubspaces);
264: for (i = 0; i < Ns; i++) PetscSpaceDestroy(&tens->tensspaces[i]);
265: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetSubspace_C", NULL);
266: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetSubspace_C", NULL);
267: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetNumSubspaces_C", NULL);
268: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetNumSubspaces_C", NULL);
269: PetscFree(tens->tensspaces);
270: PetscFree(tens);
271: return 0;
272: }
274: static PetscErrorCode PetscSpaceGetDimension_Tensor(PetscSpace sp, PetscInt *dim)
275: {
276: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
277: PetscInt i, Ns, d;
279: PetscSpaceSetUp(sp);
280: Ns = tens->numTensSpaces;
281: d = 1;
282: for (i = 0; i < Ns; i++) {
283: PetscInt id;
285: PetscSpaceGetDimension(tens->tensspaces[i], &id);
286: d *= id;
287: }
288: *dim = d;
289: return 0;
290: }
292: static PetscErrorCode PetscSpaceEvaluate_Tensor(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
293: {
294: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
295: DM dm = sp->dm;
296: PetscInt Nc = sp->Nc;
297: PetscInt Nv = sp->Nv;
298: PetscInt Ns;
299: PetscReal *lpoints, *sB = NULL, *sD = NULL, *sH = NULL;
300: PetscInt pdim;
302: if (!tens->setupCalled) {
303: PetscSpaceSetUp(sp);
304: PetscSpaceEvaluate(sp, npoints, points, B, D, H);
305: return 0;
306: }
307: Ns = tens->numTensSpaces;
308: PetscSpaceGetDimension(sp,&pdim);
309: DMGetWorkArray(dm, npoints*Nv, MPIU_REAL, &lpoints);
310: if (B || D || H) DMGetWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &sB);
311: if (D || H) DMGetWorkArray(dm, npoints*pdim*Nc*Nv, MPIU_REAL, &sD);
312: if (H) DMGetWorkArray(dm, npoints*pdim*Nc*Nv*Nv, MPIU_REAL, &sH);
313: if (B) {
314: for (PetscInt i = 0; i < npoints*pdim*Nc; i++) B[i] = 1.;
315: }
316: if (D) {
317: for (PetscInt i = 0; i < npoints*pdim*Nc*Nv; i++) D[i] = 1.;
318: }
319: if (H) {
320: for (PetscInt i = 0; i < npoints*pdim*Nc*Nv*Nv; i++) H[i] = 1.;
321: }
322: for (PetscInt s = 0, d = 0, vstep = 1, cstep = 1; s < Ns; s++) {
323: PetscInt sNv, sNc, spdim;
324: PetscInt vskip, cskip;
326: PetscSpaceGetNumVariables(tens->tensspaces[s], &sNv);
327: PetscSpaceGetNumComponents(tens->tensspaces[s], &sNc);
328: PetscSpaceGetDimension(tens->tensspaces[s], &spdim);
331: vskip = pdim / (vstep * spdim);
332: cskip = Nc / (cstep * sNc);
333: for (PetscInt p = 0; p < npoints; ++p) {
334: for (PetscInt i = 0; i < sNv; i++) {
335: lpoints[p * sNv + i] = points[p*Nv + d + i];
336: }
337: }
338: PetscSpaceEvaluate(tens->tensspaces[s], npoints, lpoints, sB, sD, sH);
339: if (B) {
340: for (PetscInt k = 0; k < vskip; k++) {
341: for (PetscInt si = 0; si < spdim; si++) {
342: for (PetscInt j = 0; j < vstep; j++) {
343: PetscInt i = (k * spdim + si) * vstep + j;
345: for (PetscInt l = 0; l < cskip; l++) {
346: for (PetscInt sc = 0; sc < sNc; sc++) {
347: for (PetscInt m = 0; m < cstep; m++) {
348: PetscInt c = (l * sNc + sc) * cstep + m;
350: for (PetscInt p = 0; p < npoints; p++) {
351: B[(pdim * p + i) * Nc + c] *= sB[(spdim * p + si) * sNc + sc];
352: }
353: }
354: }
355: }
356: }
357: }
358: }
359: }
360: if (D) {
361: for (PetscInt k = 0; k < vskip; k++) {
362: for (PetscInt si = 0; si < spdim; si++) {
363: for (PetscInt j = 0; j < vstep; j++) {
364: PetscInt i = (k * spdim + si) * vstep + j;
366: for (PetscInt l = 0; l < cskip; l++) {
367: for (PetscInt sc = 0; sc < sNc; sc++) {
368: for (PetscInt m = 0; m < cstep; m++) {
369: PetscInt c = (l * sNc + sc) * cstep + m;
371: for (PetscInt der = 0; der < Nv; der++) {
372: if (der >= d && der < d + sNv) {
373: for (PetscInt p = 0; p < npoints; p++) {
374: D[((pdim * p + i) * Nc + c)*Nv + der] *= sD[((spdim * p + si) * sNc + sc) * sNv + der - d];
375: }
376: } else {
377: for (PetscInt p = 0; p < npoints; p++) {
378: D[((pdim * p + i) * Nc + c)*Nv + der] *= sB[(spdim * p + si) * sNc + sc];
379: }
380: }
381: }
382: }
383: }
384: }
385: }
386: }
387: }
388: }
389: if (H) {
390: for (PetscInt k = 0; k < vskip; k++) {
391: for (PetscInt si = 0; si < spdim; si++) {
392: for (PetscInt j = 0; j < vstep; j++) {
393: PetscInt i = (k * spdim + si) * vstep + j;
395: for (PetscInt l = 0; l < cskip; l++) {
396: for (PetscInt sc = 0; sc < sNc; sc++) {
397: for (PetscInt m = 0; m < cstep; m++) {
398: PetscInt c = (l * sNc + sc) * cstep + m;
400: for (PetscInt der = 0; der < Nv; der++) {
401: for (PetscInt der2 = 0; der2 < Nv; der2++) {
402: if (der >= d && der < d + sNv && der2 >= d && der2 < d + sNv) {
403: for (PetscInt p = 0; p < npoints; p++) {
404: H[(((pdim * p + i) * Nc + c)*Nv + der) * Nv + der2] *= sH[(((spdim * p + si) * sNc + sc) * sNv + der - d) * sNv + der2 - d];
405: }
406: } else if (der >= d && der < d + sNv) {
407: for (PetscInt p = 0; p < npoints; p++) {
408: H[(((pdim * p + i) * Nc + c)*Nv + der) * Nv + der2] *= sD[((spdim * p + si) * sNc + sc) * sNv + der - d];
409: }
410: } else if (der2 >= d && der2 < d + sNv) {
411: for (PetscInt p = 0; p < npoints; p++) {
412: H[(((pdim * p + i) * Nc + c)*Nv + der) * Nv + der2] *= sD[((spdim * p + si) * sNc + sc) * sNv + der2 - d];
413: }
414: } else {
415: for (PetscInt p = 0; p < npoints; p++) {
416: H[(((pdim * p + i) * Nc + c)*Nv + der) * Nv + der2] *= sB[(spdim * p + si) * sNc + sc];
417: }
418: }
419: }
420: }
421: }
422: }
423: }
424: }
425: }
426: }
427: }
428: d += sNv;
429: vstep *= spdim;
430: cstep *= sNc;
431: }
432: if (H) DMRestoreWorkArray(dm, npoints*pdim*Nc*Nv*Nv, MPIU_REAL, &sH);
433: if (D || H) DMRestoreWorkArray(dm, npoints*pdim*Nc*Nv, MPIU_REAL, &sD);
434: if (B || D || H) DMRestoreWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &sB);
435: DMRestoreWorkArray(dm, npoints*Nv, MPIU_REAL, &lpoints);
436: return 0;
437: }
439: /*@
440: PetscSpaceTensorSetNumSubspaces - Set the number of spaces in the tensor product
442: Input Parameters:
443: + sp - the function space object
444: - numTensSpaces - the number of spaces
446: Level: intermediate
448: .seealso: PetscSpaceTensorGetNumSubspaces(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
449: @*/
450: PetscErrorCode PetscSpaceTensorSetNumSubspaces(PetscSpace sp, PetscInt numTensSpaces)
451: {
453: PetscTryMethod(sp,"PetscSpaceTensorSetNumSubspaces_C",(PetscSpace,PetscInt),(sp,numTensSpaces));
454: return 0;
455: }
457: /*@
458: PetscSpaceTensorGetNumSubspaces - Get the number of spaces in the tensor product
460: Input Parameter:
461: . sp - the function space object
463: Output Parameter:
464: . numTensSpaces - the number of spaces
466: Level: intermediate
468: .seealso: PetscSpaceTensorSetNumSubspaces(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
469: @*/
470: PetscErrorCode PetscSpaceTensorGetNumSubspaces(PetscSpace sp, PetscInt *numTensSpaces)
471: {
474: PetscTryMethod(sp,"PetscSpaceTensorGetNumSubspaces_C",(PetscSpace,PetscInt*),(sp,numTensSpaces));
475: return 0;
476: }
478: /*@
479: PetscSpaceTensorSetSubspace - Set a space in the tensor product
481: Input Parameters:
482: + sp - the function space object
483: . s - The space number
484: - subsp - the number of spaces
486: Level: intermediate
488: .seealso: PetscSpaceTensorGetSubspace(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
489: @*/
490: PetscErrorCode PetscSpaceTensorSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
491: {
494: PetscTryMethod(sp,"PetscSpaceTensorSetSubspace_C",(PetscSpace,PetscInt,PetscSpace),(sp,s,subsp));
495: return 0;
496: }
498: /*@
499: PetscSpaceTensorGetSubspace - Get a space in the tensor product
501: Input Parameters:
502: + sp - the function space object
503: - s - The space number
505: Output Parameter:
506: . subsp - the PetscSpace
508: Level: intermediate
510: .seealso: PetscSpaceTensorSetSubspace(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
511: @*/
512: PetscErrorCode PetscSpaceTensorGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
513: {
516: PetscTryMethod(sp,"PetscSpaceTensorGetSubspace_C",(PetscSpace,PetscInt,PetscSpace*),(sp,s,subsp));
517: return 0;
518: }
520: static PetscErrorCode PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space, PetscInt numTensSpaces)
521: {
522: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;
523: PetscInt Ns;
526: Ns = tens->numTensSpaces;
527: if (numTensSpaces == Ns) return 0;
528: if (Ns >= 0) {
529: PetscInt s;
531: for (s = 0; s < Ns; s++) PetscSpaceDestroy(&tens->tensspaces[s]);
532: PetscFree(tens->tensspaces);
533: }
534: Ns = tens->numTensSpaces = numTensSpaces;
535: PetscCalloc1(Ns, &tens->tensspaces);
536: return 0;
537: }
539: static PetscErrorCode PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space, PetscInt *numTensSpaces)
540: {
541: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;
543: *numTensSpaces = tens->numTensSpaces;
544: return 0;
545: }
547: static PetscErrorCode PetscSpaceTensorSetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace subspace)
548: {
549: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;
550: PetscInt Ns;
553: Ns = tens->numTensSpaces;
556: PetscObjectReference((PetscObject)subspace);
557: PetscSpaceDestroy(&tens->tensspaces[s]);
558: tens->tensspaces[s] = subspace;
559: return 0;
560: }
562: static PetscErrorCode PetscSpaceGetHeightSubspace_Tensor(PetscSpace sp, PetscInt height, PetscSpace *subsp)
563: {
564: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
565: PetscInt Nc, dim, order, i;
566: PetscSpace bsp;
568: PetscSpaceGetNumVariables(sp, &dim);
570: PetscSpaceGetNumComponents(sp, &Nc);
571: PetscSpaceGetDegree(sp, &order, NULL);
573: if (!tens->heightsubspaces) PetscCalloc1(dim, &tens->heightsubspaces);
574: if (height <= dim) {
575: if (!tens->heightsubspaces[height-1]) {
576: PetscSpace sub;
577: const char *name;
579: PetscSpaceTensorGetSubspace(sp, 0, &bsp);
580: PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);
581: PetscObjectGetName((PetscObject) sp, &name);
582: PetscObjectSetName((PetscObject) sub, name);
583: PetscSpaceSetType(sub, PETSCSPACETENSOR);
584: PetscSpaceSetNumComponents(sub, Nc);
585: PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);
586: PetscSpaceSetNumVariables(sub, dim-height);
587: PetscSpaceTensorSetNumSubspaces(sub, dim-height);
588: for (i = 0; i < dim - height; i++) {
589: PetscSpaceTensorSetSubspace(sub, i, bsp);
590: }
591: PetscSpaceSetUp(sub);
592: tens->heightsubspaces[height-1] = sub;
593: }
594: *subsp = tens->heightsubspaces[height-1];
595: } else {
596: *subsp = NULL;
597: }
598: return 0;
599: }
601: static PetscErrorCode PetscSpaceTensorGetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace *subspace)
602: {
603: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;
604: PetscInt Ns;
606: Ns = tens->numTensSpaces;
609: *subspace = tens->tensspaces[s];
610: return 0;
611: }
613: static PetscErrorCode PetscSpaceInitialize_Tensor(PetscSpace sp)
614: {
615: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Tensor;
616: sp->ops->setup = PetscSpaceSetUp_Tensor;
617: sp->ops->view = PetscSpaceView_Tensor;
618: sp->ops->destroy = PetscSpaceDestroy_Tensor;
619: sp->ops->getdimension = PetscSpaceGetDimension_Tensor;
620: sp->ops->evaluate = PetscSpaceEvaluate_Tensor;
621: sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Tensor;
622: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetNumSubspaces_C", PetscSpaceTensorGetNumSubspaces_Tensor);
623: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetNumSubspaces_C", PetscSpaceTensorSetNumSubspaces_Tensor);
624: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetSubspace_C", PetscSpaceTensorGetSubspace_Tensor);
625: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetSubspace_C", PetscSpaceTensorSetSubspace_Tensor);
626: return 0;
627: }
629: /*MC
630: PETSCSPACETENSOR = "tensor" - A PetscSpace object that encapsulates a tensor product space.
631: A tensor product is created of the components of the subspaces as well.
633: Level: intermediate
635: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
636: M*/
638: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Tensor(PetscSpace sp)
639: {
640: PetscSpace_Tensor *tens;
643: PetscNewLog(sp,&tens);
644: sp->data = tens;
646: tens->numTensSpaces = PETSC_DEFAULT;
648: PetscSpaceInitialize_Tensor(sp);
649: return 0;
650: }