Actual source code: spacetensor.c
petsc-3.14.6 2021-03-30
1: #include <petsc/private/petscfeimpl.h>
3: static PetscErrorCode PetscSpaceTensorCreateSubspace(PetscSpace space, PetscInt Nvs, PetscSpace *subspace)
4: {
5: PetscInt degree;
6: const char *prefix;
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, 1);
16: PetscSpaceSetDegree(*subspace, degree, PETSC_DETERMINE);
17: PetscObjectSetOptionsPrefix((PetscObject)*subspace, prefix);
18: PetscObjectAppendOptionsPrefix((PetscObject)*subspace, "subspace_");
19: return(0);
20: }
22: static PetscErrorCode PetscSpaceSetFromOptions_Tensor(PetscOptionItems *PetscOptionsObject,PetscSpace sp)
23: {
24: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
25: PetscInt Ns, Nc, i, Nv, deg;
26: PetscBool uniform = PETSC_TRUE;
27: PetscErrorCode ierr;
30: PetscSpaceGetNumVariables(sp, &Nv);
31: if (!Nv) return(0);
32: PetscSpaceGetNumComponents(sp, &Nc);
33: PetscSpaceTensorGetNumSubspaces(sp, &Ns);
34: PetscSpaceGetDegree(sp, °, NULL);
35: if (Ns > 1) {
36: PetscSpace s0;
38: PetscSpaceTensorGetSubspace(sp, 0, &s0);
39: for (i = 1; i < Ns; i++) {
40: PetscSpace si;
42: PetscSpaceTensorGetSubspace(sp, i, &si);
43: if (si != s0) {uniform = PETSC_FALSE; break;}
44: }
45: }
46: Ns = (Ns == PETSC_DEFAULT) ? PetscMax(Nv,1) : Ns;
47: PetscOptionsHead(PetscOptionsObject,"PetscSpace tensor options");
48: PetscOptionsBoundedInt("-petscspace_tensor_spaces", "The number of subspaces", "PetscSpaceTensorSetNumSubspaces", Ns, &Ns, NULL,0);
49: PetscOptionsBool("-petscspace_tensor_uniform", "Subspaces are identical", "PetscSpaceTensorSetFromOptions", uniform, &uniform, NULL);
50: PetscOptionsTail();
51: if (Ns < 0 || (Nv > 0 && Ns == 0)) SETERRQ1(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,"Cannot have a tensor space made up of %D spaces\n",Ns);
52: if (Nv > 0 && Ns > Nv) SETERRQ2(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,"Cannot have a tensor space with %D subspaces over %D variables\n", Ns, Nv);
53: if (Ns != tens->numTensSpaces) {PetscSpaceTensorSetNumSubspaces(sp, Ns);}
54: if (uniform) {
55: PetscInt Nvs = Nv / Ns;
56: PetscSpace subspace;
58: if (Nv % Ns) SETERRQ2(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_WRONG,"Cannot use %D uniform subspaces for %D variable space\n", Ns, Nv);
59: PetscSpaceTensorGetSubspace(sp, 0, &subspace);
60: if (!subspace) {PetscSpaceTensorCreateSubspace(sp, Nvs, &subspace);}
61: else {PetscObjectReference((PetscObject)subspace);}
62: PetscSpaceSetFromOptions(subspace);
63: for (i = 0; i < Ns; i++) {PetscSpaceTensorSetSubspace(sp, i, subspace);}
64: PetscSpaceDestroy(&subspace);
65: } else {
66: for (i = 0; i < Ns; i++) {
67: PetscSpace subspace;
69: PetscSpaceTensorGetSubspace(sp, i, &subspace);
70: if (!subspace) {
71: char tprefix[128];
73: PetscSpaceTensorCreateSubspace(sp, 1, &subspace);
74: PetscSNPrintf(tprefix, 128, "%d_",(int)i);
75: PetscObjectAppendOptionsPrefix((PetscObject)subspace, tprefix);
76: } else {
77: PetscObjectReference((PetscObject)subspace);
78: }
79: PetscSpaceSetFromOptions(subspace);
80: PetscSpaceTensorSetSubspace(sp, i, subspace);
81: PetscSpaceDestroy(&subspace);
82: }
83: }
84: return(0);
85: }
87: static PetscErrorCode PetscSpaceTensorView_Ascii(PetscSpace sp, PetscViewer v)
88: {
89: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
90: PetscBool uniform = PETSC_TRUE;
91: PetscInt Ns = tens->numTensSpaces, i, n;
92: PetscErrorCode ierr;
95: for (i = 1; i < Ns; i++) {
96: if (tens->tensspaces[i] != tens->tensspaces[0]) {uniform = PETSC_FALSE; break;}
97: }
98: if (uniform) {PetscViewerASCIIPrintf(v, "Tensor space of %D subspaces (all identical)\n", Ns);}
99: else {PetscViewerASCIIPrintf(v, "Tensor space of %D subspaces\n", Ns);}
100: n = uniform ? 1 : Ns;
101: for (i = 0; i < n; i++) {
102: PetscViewerASCIIPushTab(v);
103: PetscSpaceView(tens->tensspaces[i], v);
104: PetscViewerASCIIPopTab(v);
105: }
106: return(0);
107: }
109: static PetscErrorCode PetscSpaceView_Tensor(PetscSpace sp, PetscViewer viewer)
110: {
111: PetscBool iascii;
115: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
116: if (iascii) {PetscSpaceTensorView_Ascii(sp, viewer);}
117: return(0);
118: }
120: static PetscErrorCode PetscSpaceSetUp_Tensor(PetscSpace sp)
121: {
122: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
123: PetscInt Nv, Ns, i;
124: PetscBool uniform = PETSC_TRUE;
125: PetscInt deg, maxDeg;
126: PetscErrorCode ierr;
129: if (tens->setupCalled) return(0);
130: PetscSpaceGetNumVariables(sp, &Nv);
131: PetscSpaceTensorGetNumSubspaces(sp, &Ns);
132: if (Ns == PETSC_DEFAULT) {
133: Ns = Nv;
134: PetscSpaceTensorSetNumSubspaces(sp, Ns);
135: }
136: if (!Ns) {
137: if (Nv) SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zero subspaces");
138: } else {
139: PetscSpace s0;
141: if (Nv > 0 && Ns > Nv) SETERRQ2(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_OUTOFRANGE,"Cannot have a tensor space with %D subspaces over %D variables\n", Ns, Nv);
142: PetscSpaceTensorGetSubspace(sp, 0, &s0);
143: for (i = 1; i < Ns; i++) {
144: PetscSpace si;
146: PetscSpaceTensorGetSubspace(sp, i, &si);
147: if (si != s0) {uniform = PETSC_FALSE; break;}
148: }
149: if (uniform) {
150: PetscInt Nvs = Nv / Ns;
152: if (Nv % Ns) SETERRQ2(PetscObjectComm((PetscObject)sp),PETSC_ERR_ARG_WRONG,"Cannot use %D uniform subspaces for %D variable space\n", Ns, Nv);
153: if (!s0) {PetscSpaceTensorCreateSubspace(sp, Nvs, &s0);}
154: else {PetscObjectReference((PetscObject) s0);}
155: PetscSpaceSetUp(s0);
156: for (i = 0; i < Ns; i++) {PetscSpaceTensorSetSubspace(sp, i, s0);}
157: PetscSpaceDestroy(&s0);
158: } else {
159: for (i = 0 ; i < Ns; i++) {
160: PetscSpace si;
162: PetscSpaceTensorGetSubspace(sp, i, &si);
163: if (!si) {PetscSpaceTensorCreateSubspace(sp, 1, &si);}
164: else {PetscObjectReference((PetscObject) si);}
165: PetscSpaceSetUp(si);
166: PetscSpaceTensorSetSubspace(sp, i, si);
167: PetscSpaceDestroy(&si);
168: }
169: }
170: }
171: deg = PETSC_MAX_INT;
172: maxDeg = 0;
173: for (i = 0; i < Ns; i++) {
174: PetscSpace si;
175: PetscInt iDeg, iMaxDeg;
177: PetscSpaceTensorGetSubspace(sp, i, &si);
178: PetscSpaceGetDegree(si, &iDeg, &iMaxDeg);
179: deg = PetscMin(deg, iDeg);
180: maxDeg += iMaxDeg;
181: }
182: sp->degree = deg;
183: sp->maxDegree = maxDeg;
184: tens->uniform = uniform;
185: tens->setupCalled = PETSC_TRUE;
186: return(0);
187: }
189: static PetscErrorCode PetscSpaceDestroy_Tensor(PetscSpace sp)
190: {
191: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
192: PetscInt Ns, i;
193: PetscErrorCode ierr;
196: Ns = tens->numTensSpaces;
197: if (tens->heightsubspaces) {
198: PetscInt d;
200: /* sp->Nv is the spatial dimension, so it is equal to the number
201: * of subspaces on higher co-dimension points */
202: for (d = 0; d < sp->Nv; ++d) {
203: PetscSpaceDestroy(&tens->heightsubspaces[d]);
204: }
205: }
206: PetscFree(tens->heightsubspaces);
207: for (i = 0; i < Ns; i++) {PetscSpaceDestroy(&tens->tensspaces[i]);}
208: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetSubspace_C", NULL);
209: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetSubspace_C", NULL);
210: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetNumSubspaces_C", NULL);
211: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetNumSubspaces_C", NULL);
212: PetscFree(tens->tensspaces);
213: PetscFree(tens);
214: return(0);
215: }
217: static PetscErrorCode PetscSpaceGetDimension_Tensor(PetscSpace sp, PetscInt *dim)
218: {
219: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
220: PetscInt i, Ns, Nc, d;
221: PetscErrorCode ierr;
224: PetscSpaceSetUp(sp);
225: Ns = tens->numTensSpaces;
226: Nc = sp->Nc;
227: d = 1;
228: for (i = 0; i < Ns; i++) {
229: PetscInt id;
231: PetscSpaceGetDimension(tens->tensspaces[i], &id);
232: d *= id;
233: }
234: d *= Nc;
235: *dim = d;
236: return(0);
237: }
239: static PetscErrorCode PetscSpaceEvaluate_Tensor(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
240: {
241: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
242: DM dm = sp->dm;
243: PetscInt Nc = sp->Nc;
244: PetscInt Nv = sp->Nv;
245: PetscInt Ns;
246: PetscReal *lpoints, *sB = NULL, *sD = NULL, *sH = NULL;
247: PetscInt c, pdim, d, e, der, der2, i, l, si, p, s, step;
248: PetscErrorCode ierr;
251: if (!tens->setupCalled) {PetscSpaceSetUp(sp);}
252: Ns = tens->numTensSpaces;
253: PetscSpaceGetDimension(sp,&pdim);
254: pdim /= Nc;
255: DMGetWorkArray(dm, npoints*Nv, MPIU_REAL, &lpoints);
256: if (B || D || H) {DMGetWorkArray(dm, npoints*pdim, MPIU_REAL, &sB);}
257: if (D || H) {DMGetWorkArray(dm, npoints*pdim*Nv, MPIU_REAL, &sD);}
258: if (H) {DMGetWorkArray(dm, npoints*pdim*Nv*Nv, MPIU_REAL, &sH);}
259: if (B) {
260: for (i = 0; i < npoints*pdim*Nc*Nc; i++) B[i] = 0.;
261: for (i = 0; i < npoints*pdim; i++) B[i * Nc*Nc] = 1.;
262: }
263: if (D) {
264: for (i = 0; i < npoints*pdim*Nc*Nc*Nv; i++) D[i] = 0.;
265: for (i = 0; i < npoints*pdim; i++) {
266: for (l = 0; l < Nv; l++) {
267: D[i * Nc*Nc*Nv + l] = 1.;
268: }
269: }
270: }
271: if (H) {
272: for (i = 0; i < npoints*pdim*Nc*Nc*Nv*Nv; i++) H[i] = 0.;
273: for (i = 0; i < npoints*pdim; i++) {
274: for (l = 0; l < Nv*Nv; l++) {
275: H[i * Nc*Nc*Nv*Nv + l] = 1.;
276: }
277: }
278: }
279: for (s = 0, d = 0, step = 1; s < Ns; s++) {
280: PetscInt sNv, spdim;
281: PetscInt skip, j, k;
283: PetscSpaceGetNumVariables(tens->tensspaces[s], &sNv);
284: PetscSpaceGetDimension(tens->tensspaces[s], &spdim);
285: if ((pdim % step) || (pdim % spdim)) SETERRQ6(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad tensor loop: Nv %d, Ns %D, pdim %D, s %D, step %D, spdim %D", Nv, Ns, pdim, s, step, spdim);
286: skip = pdim / (step * spdim);
287: for (p = 0; p < npoints; ++p) {
288: for (i = 0; i < sNv; i++) {
289: lpoints[p * sNv + i] = points[p*Nv + d + i];
290: }
291: }
292: PetscSpaceEvaluate(tens->tensspaces[s], npoints, lpoints, sB, sD, sH);
293: if (B) {
294: for (p = 0; p < npoints; p++) {
295: for (k = 0; k < skip; k++) {
296: for (si = 0; si < spdim; si++) {
297: for (j = 0; j < step; j++) {
298: i = (k * spdim + si) * step + j;
299: B[(pdim * p + i) * Nc * Nc] *= sB[spdim * p + si];
300: }
301: }
302: }
303: }
304: }
305: if (D) {
306: for (p = 0; p < npoints; p++) {
307: for (k = 0; k < skip; k++) {
308: for (si = 0; si < spdim; si++) {
309: for (j = 0; j < step; j++) {
310: i = (k * spdim + si) * step + j;
311: for (der = 0; der < Nv; der++) {
312: if (der >= d && der < d + sNv) {
313: D[(pdim * p + i) * Nc*Nc*Nv + der] *= sD[(spdim * p + si) * sNv + der - d];
314: } else {
315: D[(pdim * p + i) * Nc*Nc*Nv + der] *= sB[spdim * p + si];
316: }
317: }
318: }
319: }
320: }
321: }
322: }
323: if (H) {
324: for (p = 0; p < npoints; p++) {
325: for (k = 0; k < skip; k++) {
326: for (si = 0; si < spdim; si++) {
327: for (j = 0; j < step; j++) {
328: i = (k * spdim + si) * step + j;
329: for (der = 0; der < Nv; der++) {
330: for (der2 = 0; der2 < Nv; der2++) {
331: if (der >= d && der < d + sNv && der2 >= d && der2 < d + sNv) {
332: H[((pdim * p + i) * Nc*Nc*Nv + der) * Nv + der2] *= sH[((spdim * p + si) * sNv + der - d) * sNv + der2 - d];
333: } else if (der >= d && der < d + sNv) {
334: H[((pdim * p + i) * Nc*Nc*Nv + der) * Nv + der2] *= sD[(spdim * p + si) * sNv + der - d];
335: } else if (der2 >= d && der2 < d + sNv) {
336: H[((pdim * p + i) * Nc*Nc*Nv + der) * Nv + der2] *= sD[(spdim * p + si) * sNv + der2 - d];
337: } else {
338: H[((pdim * p + i) * Nc*Nc*Nv + der) * Nv + der2] *= sB[spdim * p + si];
339: }
340: }
341: }
342: }
343: }
344: }
345: }
346: }
347: d += sNv;
348: step *= spdim;
349: }
350: if (B && Nc > 1) {
351: /* Make direct sum basis for multicomponent space */
352: for (p = 0; p < npoints; ++p) {
353: for (i = 0; i < pdim; ++i) {
354: for (c = 1; c < Nc; ++c) {
355: B[(p*pdim*Nc + i*Nc + c)*Nc + c] = B[(p*pdim + i)*Nc*Nc];
356: }
357: }
358: }
359: }
360: if (D && Nc > 1) {
361: /* Make direct sum basis for multicomponent space */
362: for (p = 0; p < npoints; ++p) {
363: for (i = 0; i < pdim; ++i) {
364: for (c = 1; c < Nc; ++c) {
365: for (d = 0; d < Nv; ++d) {
366: D[((p*pdim*Nc + i*Nc + c)*Nc + c)*Nv + d] = D[(p*pdim + i)*Nc*Nc*Nv + d];
367: }
368: }
369: }
370: }
371: }
372: if (H && Nc > 1) {
373: /* Make direct sum basis for multicomponent space */
374: for (p = 0; p < npoints; ++p) {
375: for (i = 0; i < pdim; ++i) {
376: for (c = 1; c < Nc; ++c) {
377: for (d = 0; d < Nv; ++d) {
378: for (e = 0; e < Nv; ++e) {
379: H[(((p*pdim*Nc + i*Nc + c)*Nc + c)*Nv + d)*Nv + e] = H[((p*pdim + i)*Nc*Nc*Nv + d)*Nv + e];
380: }
381: }
382: }
383: }
384: }
385: }
386: if (H) {DMRestoreWorkArray(dm, npoints*pdim*Nv*Nv, MPIU_REAL, &sH);}
387: if (D || H) {DMRestoreWorkArray(dm, npoints*pdim*Nv, MPIU_REAL, &sD);}
388: if (B || D || H) {DMRestoreWorkArray(dm, npoints*pdim, MPIU_REAL, &sB);}
389: DMRestoreWorkArray(dm, npoints*Nv, MPIU_REAL, &lpoints);
390: return(0);
391: }
393: /*@
394: PetscSpaceTensorSetNumSubspaces - Set the number of spaces in the tensor product
396: Input Parameters:
397: + sp - the function space object
398: - numTensSpaces - the number of spaces
400: Level: intermediate
402: .seealso: PetscSpaceTensorGetNumSubspaces(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
403: @*/
404: PetscErrorCode PetscSpaceTensorSetNumSubspaces(PetscSpace sp, PetscInt numTensSpaces)
405: {
410: PetscTryMethod(sp,"PetscSpaceTensorSetNumSubspaces_C",(PetscSpace,PetscInt),(sp,numTensSpaces));
411: return(0);
412: }
414: /*@
415: PetscSpaceTensorGetNumSubspaces - Get the number of spaces in the tensor product
417: Input Parameter:
418: . sp - the function space object
420: Output Parameter:
421: . numTensSpaces - the number of spaces
423: Level: intermediate
425: .seealso: PetscSpaceTensorSetNumSubspaces(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
426: @*/
427: PetscErrorCode PetscSpaceTensorGetNumSubspaces(PetscSpace sp, PetscInt *numTensSpaces)
428: {
434: PetscTryMethod(sp,"PetscSpaceTensorGetNumSubspaces_C",(PetscSpace,PetscInt*),(sp,numTensSpaces));
435: return(0);
436: }
438: /*@
439: PetscSpaceTensorSetSubspace - Set a space in the tensor product
441: Input Parameters:
442: + sp - the function space object
443: . s - The space number
444: - subsp - the number of spaces
446: Level: intermediate
448: .seealso: PetscSpaceTensorGetSubspace(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
449: @*/
450: PetscErrorCode PetscSpaceTensorSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
451: {
457: PetscTryMethod(sp,"PetscSpaceTensorSetSubspace_C",(PetscSpace,PetscInt,PetscSpace),(sp,s,subsp));
458: return(0);
459: }
461: /*@
462: PetscSpaceTensorGetSubspace - Get a space in the tensor product
464: Input Parameters:
465: + sp - the function space object
466: - s - The space number
468: Output Parameter:
469: . subsp - the PetscSpace
471: Level: intermediate
473: .seealso: PetscSpaceTensorSetSubspace(), PetscSpaceSetDegree(), PetscSpaceSetNumVariables()
474: @*/
475: PetscErrorCode PetscSpaceTensorGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
476: {
482: PetscTryMethod(sp,"PetscSpaceTensorGetSubspace_C",(PetscSpace,PetscInt,PetscSpace*),(sp,s,subsp));
483: return(0);
484: }
486: static PetscErrorCode PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space, PetscInt numTensSpaces)
487: {
488: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;
489: PetscInt Ns;
490: PetscErrorCode ierr;
493: if (tens->setupCalled) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Cannot change number of subspaces after setup called\n");
494: Ns = tens->numTensSpaces;
495: if (numTensSpaces == Ns) return(0);
496: if (Ns >= 0) {
497: PetscInt s;
499: for (s = 0; s < Ns; s++) {PetscSpaceDestroy(&tens->tensspaces[s]);}
500: PetscFree(tens->tensspaces);
501: }
502: Ns = tens->numTensSpaces = numTensSpaces;
503: PetscCalloc1(Ns, &tens->tensspaces);
504: return(0);
505: }
507: static PetscErrorCode PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space, PetscInt *numTensSpaces)
508: {
509: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;
512: *numTensSpaces = tens->numTensSpaces;
513: return(0);
514: }
516: static PetscErrorCode PetscSpaceTensorSetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace subspace)
517: {
518: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;
519: PetscInt Ns;
520: PetscErrorCode ierr;
523: if (tens->setupCalled) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Cannot change subspace after setup called\n");
524: Ns = tens->numTensSpaces;
525: if (Ns < 0) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSpaceTensorSetNumSubspaces() first\n");
526: if (s < 0 || s >= Ns) SETERRQ1(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_OUTOFRANGE,"Invalid subspace number %D\n",subspace);
527: PetscObjectReference((PetscObject)subspace);
528: PetscSpaceDestroy(&tens->tensspaces[s]);
529: tens->tensspaces[s] = subspace;
530: return(0);
531: }
533: static PetscErrorCode PetscSpaceGetHeightSubspace_Tensor(PetscSpace sp, PetscInt height, PetscSpace *subsp)
534: {
535: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) sp->data;
536: PetscInt Nc, dim, order, i;
537: PetscSpace bsp;
541: if (!tens->uniform) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Can only get a generic height subspace of a uniform tensor space: this tensor space is not uniform.\n");
542: PetscSpaceGetNumComponents(sp, &Nc);
543: PetscSpaceGetNumVariables(sp, &dim);
544: PetscSpaceGetDegree(sp, &order, NULL);
545: if (height > dim || height < 0) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %D for dimension %D space", height, dim);
546: if (!tens->heightsubspaces) {PetscCalloc1(dim, &tens->heightsubspaces);}
547: if (height <= dim) {
548: if (!tens->heightsubspaces[height-1]) {
549: PetscSpace sub;
550: const char *name;
552: PetscSpaceTensorGetSubspace(sp, 0, &bsp);
553: PetscSpaceCreate(PetscObjectComm((PetscObject) sp), &sub);
554: PetscObjectGetName((PetscObject) sp, &name);
555: PetscObjectSetName((PetscObject) sub, name);
556: PetscSpaceSetType(sub, PETSCSPACETENSOR);
557: PetscSpaceSetNumComponents(sub, Nc);
558: PetscSpaceSetDegree(sub, order, PETSC_DETERMINE);
559: PetscSpaceSetNumVariables(sub, dim-height);
560: PetscSpaceTensorSetNumSubspaces(sub, dim-height);
561: for (i = 0; i < dim - height; i++) {
562: PetscSpaceTensorSetSubspace(sub, i, bsp);
563: }
564: PetscSpaceSetUp(sub);
565: tens->heightsubspaces[height-1] = sub;
566: }
567: *subsp = tens->heightsubspaces[height-1];
568: } else {
569: *subsp = NULL;
570: }
571: return(0);
572: }
574: static PetscErrorCode PetscSpaceTensorGetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace *subspace)
575: {
576: PetscSpace_Tensor *tens = (PetscSpace_Tensor *) space->data;
577: PetscInt Ns;
580: Ns = tens->numTensSpaces;
581: if (Ns < 0) SETERRQ(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_WRONGSTATE,"Must call PetscSpaceTensorSetNumSubspaces() first\n");
582: if (s < 0 || s >= Ns) SETERRQ1(PetscObjectComm((PetscObject)space),PETSC_ERR_ARG_OUTOFRANGE,"Invalid subspace number %D\n",subspace);
583: *subspace = tens->tensspaces[s];
584: return(0);
585: }
587: static PetscErrorCode PetscSpaceInitialize_Tensor(PetscSpace sp)
588: {
592: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Tensor;
593: sp->ops->setup = PetscSpaceSetUp_Tensor;
594: sp->ops->view = PetscSpaceView_Tensor;
595: sp->ops->destroy = PetscSpaceDestroy_Tensor;
596: sp->ops->getdimension = PetscSpaceGetDimension_Tensor;
597: sp->ops->evaluate = PetscSpaceEvaluate_Tensor;
598: sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Tensor;
599: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetNumSubspaces_C", PetscSpaceTensorGetNumSubspaces_Tensor);
600: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetNumSubspaces_C", PetscSpaceTensorSetNumSubspaces_Tensor);
601: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorGetSubspace_C", PetscSpaceTensorGetSubspace_Tensor);
602: PetscObjectComposeFunction((PetscObject) sp, "PetscSpaceTensorSetSubspace_C", PetscSpaceTensorSetSubspace_Tensor);
603: return(0);
604: }
606: /*MC
607: PETSCSPACETENSOR = "tensor" - A PetscSpace object that encapsulates a tensor product space.
608: Subspaces are scalar spaces (num of componenents = 1), so the components
609: of a vector-valued tensor space are assumed to be identical.
611: Level: intermediate
613: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
614: M*/
616: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Tensor(PetscSpace sp)
617: {
618: PetscSpace_Tensor *tens;
619: PetscErrorCode ierr;
623: PetscNewLog(sp,&tens);
624: sp->data = tens;
626: tens->numTensSpaces = PETSC_DEFAULT;
628: PetscSpaceInitialize_Tensor(sp);
629: return(0);
630: }