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 PetscObjectSetName((PetscObject)*subspace, "tensor component");
24: return 0;
25: }
27: static PetscErrorCode PetscSpaceSetFromOptions_Tensor(PetscSpace sp, PetscOptionItems *PetscOptionsObject)
28: {
29: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
30: PetscInt Ns, Nc, i, Nv, deg;
31: PetscBool uniform = PETSC_TRUE;
33: PetscSpaceGetNumVariables(sp, &Nv);
34: if (!Nv) return 0;
35: PetscSpaceGetNumComponents(sp, &Nc);
36: PetscSpaceTensorGetNumSubspaces(sp, &Ns);
37: PetscSpaceGetDegree(sp, °, NULL);
38: if (Ns > 1) {
39: PetscSpace s0;
41: PetscSpaceTensorGetSubspace(sp, 0, &s0);
42: for (i = 1; i < Ns; i++) {
43: PetscSpace si;
45: PetscSpaceTensorGetSubspace(sp, i, &si);
46: if (si != s0) {
47: uniform = PETSC_FALSE;
48: break;
49: }
50: }
51: }
52: Ns = (Ns == PETSC_DEFAULT) ? PetscMax(Nv, 1) : Ns;
53: PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace tensor options");
54: PetscOptionsBoundedInt("-petscspace_tensor_spaces", "The number of subspaces", "PetscSpaceTensorSetNumSubspaces", Ns, &Ns, NULL, 0);
55: PetscOptionsBool("-petscspace_tensor_uniform", "Subspaces are identical", "PetscSpaceTensorSetFromOptions", uniform, &uniform, NULL);
56: PetscOptionsHeadEnd();
59: if (Ns != tens->numTensSpaces) PetscSpaceTensorSetNumSubspaces(sp, Ns);
60: if (uniform) {
61: PetscInt Nvs = Nv / Ns;
62: PetscInt Ncs;
63: PetscSpace subspace;
66: Ncs = (PetscInt)PetscPowReal((PetscReal)Nc, 1. / Ns);
68: PetscSpaceTensorGetSubspace(sp, 0, &subspace);
69: if (!subspace) PetscSpaceTensorCreateSubspace(sp, Nvs, Ncs, &subspace);
70: else PetscObjectReference((PetscObject)subspace);
71: PetscSpaceSetFromOptions(subspace);
72: for (i = 0; i < Ns; i++) PetscSpaceTensorSetSubspace(sp, i, subspace);
73: PetscSpaceDestroy(&subspace);
74: } else {
75: for (i = 0; i < Ns; i++) {
76: PetscSpace subspace;
78: PetscSpaceTensorGetSubspace(sp, i, &subspace);
79: if (!subspace) {
80: char tprefix[128];
82: PetscSpaceTensorCreateSubspace(sp, 1, 1, &subspace);
83: PetscSNPrintf(tprefix, 128, "%d_", (int)i);
84: PetscObjectAppendOptionsPrefix((PetscObject)subspace, tprefix);
85: } else PetscObjectReference((PetscObject)subspace);
86: PetscSpaceSetFromOptions(subspace);
87: PetscSpaceTensorSetSubspace(sp, i, subspace);
88: PetscSpaceDestroy(&subspace);
89: }
90: }
91: return 0;
92: }
94: static PetscErrorCode PetscSpaceTensorView_Ascii(PetscSpace sp, PetscViewer v)
95: {
96: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
97: PetscBool uniform = PETSC_TRUE;
98: PetscInt Ns = tens->numTensSpaces, i, n;
100: for (i = 1; i < Ns; i++) {
101: if (tens->tensspaces[i] != tens->tensspaces[0]) {
102: uniform = PETSC_FALSE;
103: break;
104: }
105: }
106: if (uniform) PetscViewerASCIIPrintf(v, "Tensor space of %" PetscInt_FMT " subspaces (all identical)\n", Ns);
107: else PetscViewerASCIIPrintf(v, "Tensor space of %" PetscInt_FMT " subspaces\n", Ns);
108: n = uniform ? 1 : Ns;
109: for (i = 0; i < n; i++) {
110: PetscViewerASCIIPushTab(v);
111: PetscSpaceView(tens->tensspaces[i], v);
112: PetscViewerASCIIPopTab(v);
113: }
114: return 0;
115: }
117: static PetscErrorCode PetscSpaceView_Tensor(PetscSpace sp, PetscViewer viewer)
118: {
119: PetscBool iascii;
121: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
122: if (iascii) PetscSpaceTensorView_Ascii(sp, viewer);
123: return 0;
124: }
126: static PetscErrorCode PetscSpaceSetUp_Tensor(PetscSpace sp)
127: {
128: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
129: PetscInt Nc, Nv, Ns;
130: PetscBool uniform = PETSC_TRUE;
131: PetscInt deg, maxDeg;
132: PetscInt Ncprod;
134: if (tens->setupCalled) return 0;
135: PetscSpaceGetNumVariables(sp, &Nv);
136: PetscSpaceGetNumComponents(sp, &Nc);
137: PetscSpaceTensorGetNumSubspaces(sp, &Ns);
138: if (Ns == PETSC_DEFAULT) {
139: Ns = Nv;
140: PetscSpaceTensorSetNumSubspaces(sp, Ns);
141: }
142: if (!Ns) {
143: SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zero subspaces");
144: } else {
145: PetscSpace s0;
148: PetscSpaceTensorGetSubspace(sp, 0, &s0);
149: for (PetscInt i = 1; i < Ns; i++) {
150: PetscSpace si;
152: PetscSpaceTensorGetSubspace(sp, i, &si);
153: if (si != s0) {
154: uniform = PETSC_FALSE;
155: break;
156: }
157: }
158: if (uniform) {
159: PetscInt Nvs = Nv / Ns;
160: PetscInt Ncs;
163: Ncs = (PetscInt)(PetscPowReal((PetscReal)Nc, 1. / Ns));
165: if (!s0) PetscSpaceTensorCreateSubspace(sp, Nvs, Ncs, &s0);
166: else PetscObjectReference((PetscObject)s0);
167: PetscSpaceSetUp(s0);
168: for (PetscInt i = 0; i < Ns; i++) PetscSpaceTensorSetSubspace(sp, i, s0);
169: PetscSpaceDestroy(&s0);
170: Ncprod = PetscPowInt(Ncs, Ns);
171: } else {
172: PetscInt Nvsum = 0;
174: Ncprod = 1;
175: for (PetscInt i = 0; i < Ns; i++) {
176: PetscInt Nvs, Ncs;
177: PetscSpace si;
179: PetscSpaceTensorGetSubspace(sp, i, &si);
180: if (!si) PetscSpaceTensorCreateSubspace(sp, 1, 1, &si);
181: else PetscObjectReference((PetscObject)si);
182: PetscSpaceSetUp(si);
183: PetscSpaceTensorSetSubspace(sp, i, si);
184: PetscSpaceGetNumVariables(si, &Nvs);
185: PetscSpaceGetNumComponents(si, &Ncs);
186: PetscSpaceDestroy(&si);
188: Nvsum += Nvs;
189: Ncprod *= Ncs;
190: }
194: }
195: if (Ncprod != Nc) {
196: PetscInt Ncopies = Nc / Ncprod;
197: PetscInt Nv = sp->Nv;
198: const char *prefix;
199: const char *name;
200: char subname[PETSC_MAX_PATH_LEN];
201: PetscSpace subsp;
203: PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp);
204: PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix);
205: PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix);
206: PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_");
207: if (((PetscObject)sp)->name) {
208: PetscObjectGetName((PetscObject)sp, &name);
209: PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s sum component", name);
210: PetscObjectSetName((PetscObject)subsp, subname);
211: } else PetscObjectSetName((PetscObject)subsp, "sum component");
212: PetscSpaceSetType(subsp, PETSCSPACETENSOR);
213: PetscSpaceSetNumVariables(subsp, Nv);
214: PetscSpaceSetNumComponents(subsp, Ncprod);
215: PetscSpaceTensorSetNumSubspaces(subsp, Ns);
216: for (PetscInt i = 0; i < Ns; i++) {
217: PetscSpace ssp;
219: PetscSpaceTensorGetSubspace(sp, i, &ssp);
220: PetscSpaceTensorSetSubspace(subsp, i, ssp);
221: }
222: PetscSpaceSetUp(subsp);
223: PetscSpaceSetType(sp, PETSCSPACESUM);
224: PetscSpaceSumSetNumSubspaces(sp, Ncopies);
225: for (PetscInt i = 0; i < Ncopies; i++) PetscSpaceSumSetSubspace(sp, i, subsp);
226: PetscSpaceDestroy(&subsp);
227: PetscSpaceSetUp(sp);
228: return 0;
229: }
230: }
231: deg = PETSC_MAX_INT;
232: maxDeg = 0;
233: for (PetscInt i = 0; i < Ns; i++) {
234: PetscSpace si;
235: PetscInt iDeg, iMaxDeg;
237: PetscSpaceTensorGetSubspace(sp, i, &si);
238: PetscSpaceGetDegree(si, &iDeg, &iMaxDeg);
239: deg = PetscMin(deg, iDeg);
240: maxDeg += iMaxDeg;
241: }
242: sp->degree = deg;
243: sp->maxDegree = maxDeg;
244: tens->uniform = uniform;
245: tens->setupCalled = PETSC_TRUE;
246: return 0;
247: }
249: static PetscErrorCode PetscSpaceDestroy_Tensor(PetscSpace sp)
250: {
251: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
252: PetscInt Ns, i;
254: Ns = tens->numTensSpaces;
255: if (tens->heightsubspaces) {
256: PetscInt d;
258: /* sp->Nv is the spatial dimension, so it is equal to the number
259: * of subspaces on higher co-dimension points */
260: for (d = 0; d < sp->Nv; ++d) PetscSpaceDestroy(&tens->heightsubspaces[d]);
261: }
262: PetscFree(tens->heightsubspaces);
263: for (i = 0; i < Ns; i++) PetscSpaceDestroy(&tens->tensspaces[i]);
264: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", NULL);
265: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", NULL);
266: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", NULL);
267: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", NULL);
268: PetscFree(tens->tensspaces);
269: PetscFree(tens);
270: return 0;
271: }
273: static PetscErrorCode PetscSpaceGetDimension_Tensor(PetscSpace sp, PetscInt *dim)
274: {
275: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
276: PetscInt i, Ns, d;
278: PetscSpaceSetUp(sp);
279: Ns = tens->numTensSpaces;
280: d = 1;
281: for (i = 0; i < Ns; i++) {
282: PetscInt id;
284: PetscSpaceGetDimension(tens->tensspaces[i], &id);
285: d *= id;
286: }
287: *dim = d;
288: return 0;
289: }
291: static PetscErrorCode PetscSpaceEvaluate_Tensor(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
292: {
293: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
294: DM dm = sp->dm;
295: PetscInt Nc = sp->Nc;
296: PetscInt Nv = sp->Nv;
297: PetscInt Ns;
298: PetscReal *lpoints, *sB = NULL, *sD = NULL, *sH = NULL;
299: PetscInt pdim;
301: if (!tens->setupCalled) {
302: PetscSpaceSetUp(sp);
303: PetscSpaceEvaluate(sp, npoints, points, B, D, H);
304: return 0;
305: }
306: Ns = tens->numTensSpaces;
307: PetscSpaceGetDimension(sp, &pdim);
308: DMGetWorkArray(dm, npoints * Nv, MPIU_REAL, &lpoints);
309: if (B || D || H) DMGetWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &sB);
310: if (D || H) DMGetWorkArray(dm, npoints * pdim * Nc * Nv, MPIU_REAL, &sD);
311: if (H) DMGetWorkArray(dm, npoints * pdim * Nc * Nv * Nv, MPIU_REAL, &sH);
312: if (B) {
313: for (PetscInt i = 0; i < npoints * pdim * Nc; i++) B[i] = 1.;
314: }
315: if (D) {
316: for (PetscInt i = 0; i < npoints * pdim * Nc * Nv; i++) D[i] = 1.;
317: }
318: if (H) {
319: for (PetscInt i = 0; i < npoints * pdim * Nc * Nv * Nv; i++) H[i] = 1.;
320: }
321: for (PetscInt s = 0, d = 0, vstep = 1, cstep = 1; s < Ns; s++) {
322: PetscInt sNv, sNc, spdim;
323: PetscInt vskip, cskip;
325: PetscSpaceGetNumVariables(tens->tensspaces[s], &sNv);
326: PetscSpaceGetNumComponents(tens->tensspaces[s], &sNc);
327: PetscSpaceGetDimension(tens->tensspaces[s], &spdim);
330: vskip = pdim / (vstep * spdim);
331: cskip = Nc / (cstep * sNc);
332: for (PetscInt p = 0; p < npoints; ++p) {
333: for (PetscInt i = 0; i < sNv; i++) lpoints[p * sNv + i] = points[p * Nv + d + i];
334: }
335: PetscSpaceEvaluate(tens->tensspaces[s], npoints, lpoints, sB, sD, sH);
336: if (B) {
337: for (PetscInt k = 0; k < vskip; k++) {
338: for (PetscInt si = 0; si < spdim; si++) {
339: for (PetscInt j = 0; j < vstep; j++) {
340: PetscInt i = (k * spdim + si) * vstep + j;
342: for (PetscInt l = 0; l < cskip; l++) {
343: for (PetscInt sc = 0; sc < sNc; sc++) {
344: for (PetscInt m = 0; m < cstep; m++) {
345: PetscInt c = (l * sNc + sc) * cstep + m;
347: for (PetscInt p = 0; p < npoints; p++) B[(pdim * p + i) * Nc + c] *= sB[(spdim * p + si) * sNc + sc];
348: }
349: }
350: }
351: }
352: }
353: }
354: }
355: if (D) {
356: for (PetscInt k = 0; k < vskip; k++) {
357: for (PetscInt si = 0; si < spdim; si++) {
358: for (PetscInt j = 0; j < vstep; j++) {
359: PetscInt i = (k * spdim + si) * vstep + j;
361: for (PetscInt l = 0; l < cskip; l++) {
362: for (PetscInt sc = 0; sc < sNc; sc++) {
363: for (PetscInt m = 0; m < cstep; m++) {
364: PetscInt c = (l * sNc + sc) * cstep + m;
366: for (PetscInt der = 0; der < Nv; der++) {
367: if (der >= d && der < d + sNv) {
368: for (PetscInt p = 0; p < npoints; p++) D[((pdim * p + i) * Nc + c) * Nv + der] *= sD[((spdim * p + si) * sNc + sc) * sNv + der - d];
369: } else {
370: for (PetscInt p = 0; p < npoints; p++) D[((pdim * p + i) * Nc + c) * Nv + der] *= sB[(spdim * p + si) * sNc + sc];
371: }
372: }
373: }
374: }
375: }
376: }
377: }
378: }
379: }
380: if (H) {
381: for (PetscInt k = 0; k < vskip; k++) {
382: for (PetscInt si = 0; si < spdim; si++) {
383: for (PetscInt j = 0; j < vstep; j++) {
384: PetscInt i = (k * spdim + si) * vstep + j;
386: for (PetscInt l = 0; l < cskip; l++) {
387: for (PetscInt sc = 0; sc < sNc; sc++) {
388: for (PetscInt m = 0; m < cstep; m++) {
389: PetscInt c = (l * sNc + sc) * cstep + m;
391: for (PetscInt der = 0; der < Nv; der++) {
392: for (PetscInt der2 = 0; der2 < Nv; der2++) {
393: if (der >= d && der < d + sNv && der2 >= d && der2 < d + sNv) {
394: for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sH[(((spdim * p + si) * sNc + sc) * sNv + der - d) * sNv + der2 - d];
395: } else if (der >= d && der < d + sNv) {
396: for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sD[((spdim * p + si) * sNc + sc) * sNv + der - d];
397: } else if (der2 >= d && der2 < d + sNv) {
398: for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sD[((spdim * p + si) * sNc + sc) * sNv + der2 - d];
399: } else {
400: for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sB[(spdim * p + si) * sNc + sc];
401: }
402: }
403: }
404: }
405: }
406: }
407: }
408: }
409: }
410: }
411: d += sNv;
412: vstep *= spdim;
413: cstep *= sNc;
414: }
415: if (H) DMRestoreWorkArray(dm, npoints * pdim * Nc * Nv * Nv, MPIU_REAL, &sH);
416: if (D || H) DMRestoreWorkArray(dm, npoints * pdim * Nc * Nv, MPIU_REAL, &sD);
417: if (B || D || H) DMRestoreWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &sB);
418: DMRestoreWorkArray(dm, npoints * Nv, MPIU_REAL, &lpoints);
419: return 0;
420: }
422: /*@
423: PetscSpaceTensorSetNumSubspaces - Set the number of spaces in the tensor product
425: Input Parameters:
426: + sp - the function space object
427: - numTensSpaces - the number of spaces
429: Level: intermediate
431: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
432: @*/
433: PetscErrorCode PetscSpaceTensorSetNumSubspaces(PetscSpace sp, PetscInt numTensSpaces)
434: {
436: PetscTryMethod(sp, "PetscSpaceTensorSetNumSubspaces_C", (PetscSpace, PetscInt), (sp, numTensSpaces));
437: return 0;
438: }
440: /*@
441: PetscSpaceTensorGetNumSubspaces - Get the number of spaces in the tensor product
443: Input Parameter:
444: . sp - the function space object
446: Output Parameter:
447: . numTensSpaces - the number of spaces
449: Level: intermediate
451: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
452: @*/
453: PetscErrorCode PetscSpaceTensorGetNumSubspaces(PetscSpace sp, PetscInt *numTensSpaces)
454: {
457: PetscTryMethod(sp, "PetscSpaceTensorGetNumSubspaces_C", (PetscSpace, PetscInt *), (sp, numTensSpaces));
458: return 0;
459: }
461: /*@
462: PetscSpaceTensorSetSubspace - Set a space in the tensor product
464: Input Parameters:
465: + sp - the function space object
466: . s - The space number
467: - subsp - the number of spaces
469: Level: intermediate
471: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
472: @*/
473: PetscErrorCode PetscSpaceTensorSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
474: {
477: PetscTryMethod(sp, "PetscSpaceTensorSetSubspace_C", (PetscSpace, PetscInt, PetscSpace), (sp, s, subsp));
478: return 0;
479: }
481: /*@
482: PetscSpaceTensorGetSubspace - Get a space in the tensor product
484: Input Parameters:
485: + sp - the function space object
486: - s - The space number
488: Output Parameter:
489: . subsp - the PetscSpace
491: Level: intermediate
493: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
494: @*/
495: PetscErrorCode PetscSpaceTensorGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
496: {
499: PetscTryMethod(sp, "PetscSpaceTensorGetSubspace_C", (PetscSpace, PetscInt, PetscSpace *), (sp, s, subsp));
500: return 0;
501: }
503: static PetscErrorCode PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space, PetscInt numTensSpaces)
504: {
505: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
506: PetscInt Ns;
509: Ns = tens->numTensSpaces;
510: if (numTensSpaces == Ns) return 0;
511: if (Ns >= 0) {
512: PetscInt s;
514: for (s = 0; s < Ns; s++) PetscSpaceDestroy(&tens->tensspaces[s]);
515: PetscFree(tens->tensspaces);
516: }
517: Ns = tens->numTensSpaces = numTensSpaces;
518: PetscCalloc1(Ns, &tens->tensspaces);
519: return 0;
520: }
522: static PetscErrorCode PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space, PetscInt *numTensSpaces)
523: {
524: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
526: *numTensSpaces = tens->numTensSpaces;
527: return 0;
528: }
530: static PetscErrorCode PetscSpaceTensorSetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace subspace)
531: {
532: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
533: PetscInt Ns;
536: Ns = tens->numTensSpaces;
539: PetscObjectReference((PetscObject)subspace);
540: PetscSpaceDestroy(&tens->tensspaces[s]);
541: tens->tensspaces[s] = subspace;
542: return 0;
543: }
545: static PetscErrorCode PetscSpaceGetHeightSubspace_Tensor(PetscSpace sp, PetscInt height, PetscSpace *subsp)
546: {
547: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
548: PetscInt Nc, dim, order, i;
549: PetscSpace bsp;
551: PetscSpaceGetNumVariables(sp, &dim);
553: PetscSpaceGetNumComponents(sp, &Nc);
554: PetscSpaceGetDegree(sp, &order, NULL);
556: if (!tens->heightsubspaces) PetscCalloc1(dim, &tens->heightsubspaces);
557: if (height <= dim) {
558: if (!tens->heightsubspaces[height - 1]) {
559: PetscSpace sub;
560: const char *name;
562: PetscSpaceTensorGetSubspace(sp, 0, &bsp);
563: PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub);
564: PetscObjectGetName((PetscObject)sp, &name);
565: PetscObjectSetName((PetscObject)sub, name);
566: PetscSpaceSetType(sub, PETSCSPACETENSOR);
567: PetscSpaceSetNumComponents(sub, Nc);
568: PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);
569: PetscSpaceSetNumVariables(sub, dim - height);
570: PetscSpaceTensorSetNumSubspaces(sub, dim - height);
571: for (i = 0; i < dim - height; i++) PetscSpaceTensorSetSubspace(sub, i, bsp);
572: PetscSpaceSetUp(sub);
573: tens->heightsubspaces[height - 1] = sub;
574: }
575: *subsp = tens->heightsubspaces[height - 1];
576: } else {
577: *subsp = NULL;
578: }
579: return 0;
580: }
582: static PetscErrorCode PetscSpaceTensorGetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace *subspace)
583: {
584: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
585: PetscInt Ns;
587: Ns = tens->numTensSpaces;
590: *subspace = tens->tensspaces[s];
591: return 0;
592: }
594: static PetscErrorCode PetscSpaceInitialize_Tensor(PetscSpace sp)
595: {
596: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Tensor;
597: sp->ops->setup = PetscSpaceSetUp_Tensor;
598: sp->ops->view = PetscSpaceView_Tensor;
599: sp->ops->destroy = PetscSpaceDestroy_Tensor;
600: sp->ops->getdimension = PetscSpaceGetDimension_Tensor;
601: sp->ops->evaluate = PetscSpaceEvaluate_Tensor;
602: sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Tensor;
603: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", PetscSpaceTensorGetNumSubspaces_Tensor);
604: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", PetscSpaceTensorSetNumSubspaces_Tensor);
605: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", PetscSpaceTensorGetSubspace_Tensor);
606: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", PetscSpaceTensorSetSubspace_Tensor);
607: return 0;
608: }
610: /*MC
611: PETSCSPACETENSOR = "tensor" - A `PetscSpace` object that encapsulates a tensor product space.
612: A tensor product is created of the components of the subspaces as well.
614: Level: intermediate
616: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
617: M*/
619: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Tensor(PetscSpace sp)
620: {
621: PetscSpace_Tensor *tens;
624: PetscNew(&tens);
625: sp->data = tens;
627: tens->numTensSpaces = PETSC_DEFAULT;
629: PetscSpaceInitialize_Tensor(sp);
630: return 0;
631: }