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: PetscFunctionBegin;
11: PetscCall(PetscSpaceGetDegree(space, °ree, NULL));
12: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)space, &prefix));
13: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)space), subspace));
14: PetscCall(PetscSpaceSetType(*subspace, PETSCSPACEPOLYNOMIAL));
15: PetscCall(PetscSpaceSetNumVariables(*subspace, Nvs));
16: PetscCall(PetscSpaceSetNumComponents(*subspace, Ncs));
17: PetscCall(PetscSpaceSetDegree(*subspace, degree, PETSC_DETERMINE));
18: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)*subspace, prefix));
19: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)*subspace, "tensorcomp_"));
20: if (((PetscObject)space)->name) {
21: PetscCall(PetscObjectGetName((PetscObject)space, &name));
22: PetscCall(PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s tensor component", name));
23: PetscCall(PetscObjectSetName((PetscObject)*subspace, subname));
24: } else PetscCall(PetscObjectSetName((PetscObject)*subspace, "tensor component"));
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: static PetscErrorCode PetscSpaceSetFromOptions_Tensor(PetscSpace sp, PetscOptionItems PetscOptionsObject)
29: {
30: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
31: PetscInt Ns, Nc, i, Nv, deg;
32: PetscBool uniform = PETSC_TRUE;
34: PetscFunctionBegin;
35: PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
36: if (!Nv) PetscFunctionReturn(PETSC_SUCCESS);
37: PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
38: PetscCall(PetscSpaceTensorGetNumSubspaces(sp, &Ns));
39: PetscCall(PetscSpaceGetDegree(sp, °, NULL));
40: if (Ns > 1) {
41: PetscSpace s0;
43: PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &s0));
44: for (i = 1; i < Ns; i++) {
45: PetscSpace si;
47: PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si));
48: if (si != s0) {
49: uniform = PETSC_FALSE;
50: break;
51: }
52: }
53: }
54: Ns = (Ns == PETSC_DEFAULT) ? PetscMax(Nv, 1) : Ns;
55: PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace tensor options");
56: PetscCall(PetscOptionsBoundedInt("-petscspace_tensor_spaces", "The number of subspaces", "PetscSpaceTensorSetNumSubspaces", Ns, &Ns, NULL, 0));
57: PetscCall(PetscOptionsBool("-petscspace_tensor_uniform", "Subspaces are identical", "PetscSpaceTensorSetFromOptions", uniform, &uniform, NULL));
58: PetscOptionsHeadEnd();
59: PetscCheck(Ns >= 0 && (Nv <= 0 || Ns != 0), PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a tensor space made up of %" PetscInt_FMT " spaces", Ns);
60: PetscCheck(Nv <= 0 || Ns <= Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a tensor space with %" PetscInt_FMT " subspaces over %" PetscInt_FMT " variables", Ns, Nv);
61: if (Ns != tens->numTensSpaces) PetscCall(PetscSpaceTensorSetNumSubspaces(sp, Ns));
62: if (uniform) {
63: PetscInt Nvs = Nv / Ns;
64: PetscInt Ncs;
65: PetscSpace subspace;
67: PetscCheck(Nv % Ns == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " variable space", Ns, Nv);
68: Ncs = (PetscInt)PetscPowReal((PetscReal)Nc, 1. / Ns);
69: PetscCheck(Nc % PetscPowInt(Ncs, Ns) == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " component space", Ns, Nc);
70: PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &subspace));
71: if (!subspace) PetscCall(PetscSpaceTensorCreateSubspace(sp, Nvs, Ncs, &subspace));
72: else PetscCall(PetscObjectReference((PetscObject)subspace));
73: PetscCall(PetscSpaceSetFromOptions(subspace));
74: for (i = 0; i < Ns; i++) PetscCall(PetscSpaceTensorSetSubspace(sp, i, subspace));
75: PetscCall(PetscSpaceDestroy(&subspace));
76: } else {
77: for (i = 0; i < Ns; i++) {
78: PetscSpace subspace;
80: PetscCall(PetscSpaceTensorGetSubspace(sp, i, &subspace));
81: if (!subspace) {
82: char tprefix[128];
84: PetscCall(PetscSpaceTensorCreateSubspace(sp, 1, 1, &subspace));
85: PetscCall(PetscSNPrintf(tprefix, 128, "%" PetscInt_FMT "_", i));
86: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subspace, tprefix));
87: } else PetscCall(PetscObjectReference((PetscObject)subspace));
88: PetscCall(PetscSpaceSetFromOptions(subspace));
89: PetscCall(PetscSpaceTensorSetSubspace(sp, i, subspace));
90: PetscCall(PetscSpaceDestroy(&subspace));
91: }
92: }
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: static PetscErrorCode PetscSpaceTensorView_Ascii(PetscSpace sp, PetscViewer v)
97: {
98: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
99: PetscBool uniform = PETSC_TRUE;
100: PetscInt Ns = tens->numTensSpaces, i, n;
102: PetscFunctionBegin;
103: for (i = 1; i < Ns; i++) {
104: if (tens->tensspaces[i] != tens->tensspaces[0]) {
105: uniform = PETSC_FALSE;
106: break;
107: }
108: }
109: if (uniform) PetscCall(PetscViewerASCIIPrintf(v, "Tensor space of %" PetscInt_FMT " subspaces (all identical)\n", Ns));
110: else PetscCall(PetscViewerASCIIPrintf(v, "Tensor space of %" PetscInt_FMT " subspaces\n", Ns));
111: n = uniform ? 1 : Ns;
112: for (i = 0; i < n; i++) {
113: PetscCall(PetscViewerASCIIPushTab(v));
114: PetscCall(PetscSpaceView(tens->tensspaces[i], v));
115: PetscCall(PetscViewerASCIIPopTab(v));
116: }
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: static PetscErrorCode PetscSpaceView_Tensor(PetscSpace sp, PetscViewer viewer)
121: {
122: PetscBool isascii;
124: PetscFunctionBegin;
125: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
126: if (isascii) PetscCall(PetscSpaceTensorView_Ascii(sp, viewer));
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: static PetscErrorCode PetscSpaceSetUp_Tensor(PetscSpace sp)
131: {
132: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
133: PetscInt Nc, Nv, Ns;
134: PetscBool uniform = PETSC_TRUE;
135: PetscInt deg, maxDeg;
136: PetscInt Ncprod;
138: PetscFunctionBegin;
139: if (tens->setupcalled) PetscFunctionReturn(PETSC_SUCCESS);
140: PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
141: PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
142: PetscCall(PetscSpaceTensorGetNumSubspaces(sp, &Ns));
143: if (Ns == PETSC_DEFAULT) {
144: Ns = Nv;
145: PetscCall(PetscSpaceTensorSetNumSubspaces(sp, Ns));
146: }
147: if (!Ns) {
148: SETERRQ(PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have zero subspaces");
149: } else {
150: PetscSpace s0 = NULL;
152: PetscCheck(Nv <= 0 || Ns <= Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a tensor space with %" PetscInt_FMT " subspaces over %" PetscInt_FMT " variables", Ns, Nv);
153: PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &s0));
154: for (PetscInt i = 1; i < Ns; i++) {
155: PetscSpace si;
157: PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si));
158: if (si != s0) {
159: uniform = PETSC_FALSE;
160: break;
161: }
162: }
163: if (uniform) {
164: PetscInt Nvs = Nv / Ns;
165: PetscInt Ncs;
167: PetscCheck(Nv % Ns == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " variable space", Ns, Nv);
168: Ncs = (PetscInt)(PetscPowReal((PetscReal)Nc, 1. / Ns));
169: PetscCheck(Nc % PetscPowInt(Ncs, Ns) == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Cannot use %" PetscInt_FMT " uniform subspaces for %" PetscInt_FMT " component space", Ns, Nc);
170: if (!s0) PetscCall(PetscSpaceTensorCreateSubspace(sp, Nvs, Ncs, &s0));
171: else PetscCall(PetscObjectReference((PetscObject)s0));
172: PetscCall(PetscSpaceSetUp(s0));
173: for (PetscInt i = 0; i < Ns; i++) PetscCall(PetscSpaceTensorSetSubspace(sp, i, s0));
174: PetscCall(PetscSpaceDestroy(&s0));
175: Ncprod = PetscPowInt(Ncs, Ns);
176: } else {
177: PetscInt Nvsum = 0;
179: Ncprod = 1;
180: for (PetscInt i = 0; i < Ns; i++) {
181: PetscInt Nvs, Ncs;
182: PetscSpace si = NULL;
184: PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si));
185: if (!si) PetscCall(PetscSpaceTensorCreateSubspace(sp, 1, 1, &si));
186: else PetscCall(PetscObjectReference((PetscObject)si));
187: PetscCall(PetscSpaceSetUp(si));
188: PetscCall(PetscSpaceTensorSetSubspace(sp, i, si));
189: PetscCall(PetscSpaceGetNumVariables(si, &Nvs));
190: PetscCall(PetscSpaceGetNumComponents(si, &Ncs));
191: PetscCall(PetscSpaceDestroy(&si));
193: Nvsum += Nvs;
194: Ncprod *= Ncs;
195: }
197: PetscCheck(Nvsum == Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Sum of subspace variables %" PetscInt_FMT " does not equal the number of variables %" PetscInt_FMT, Nvsum, Nv);
198: PetscCheck(Nc % Ncprod == 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONG, "Product of subspace components %" PetscInt_FMT " does not divide the number of components %" PetscInt_FMT, Ncprod, Nc);
199: }
200: if (Ncprod != Nc) {
201: PetscInt Ncopies = Nc / Ncprod;
202: PetscInt Nv = sp->Nv;
203: const char *prefix;
204: const char *name;
205: char subname[PETSC_MAX_PATH_LEN];
206: PetscSpace subsp;
208: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subsp));
209: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix));
210: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subsp, prefix));
211: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subsp, "sumcomp_"));
212: if (((PetscObject)sp)->name) {
213: PetscCall(PetscObjectGetName((PetscObject)sp, &name));
214: PetscCall(PetscSNPrintf(subname, PETSC_MAX_PATH_LEN - 1, "%s sum component", name));
215: PetscCall(PetscObjectSetName((PetscObject)subsp, subname));
216: } else PetscCall(PetscObjectSetName((PetscObject)subsp, "sum component"));
217: PetscCall(PetscSpaceSetType(subsp, PETSCSPACETENSOR));
218: PetscCall(PetscSpaceSetNumVariables(subsp, Nv));
219: PetscCall(PetscSpaceSetNumComponents(subsp, Ncprod));
220: PetscCall(PetscSpaceTensorSetNumSubspaces(subsp, Ns));
221: for (PetscInt i = 0; i < Ns; i++) {
222: PetscSpace ssp;
224: PetscCall(PetscSpaceTensorGetSubspace(sp, i, &ssp));
225: PetscCall(PetscSpaceTensorSetSubspace(subsp, i, ssp));
226: }
227: PetscCall(PetscSpaceSetUp(subsp));
228: PetscCall(PetscSpaceSetType(sp, PETSCSPACESUM));
229: PetscCall(PetscSpaceSumSetNumSubspaces(sp, Ncopies));
230: for (PetscInt i = 0; i < Ncopies; i++) PetscCall(PetscSpaceSumSetSubspace(sp, i, subsp));
231: PetscCall(PetscSpaceDestroy(&subsp));
232: PetscCall(PetscSpaceSetUp(sp));
233: PetscFunctionReturn(PETSC_SUCCESS);
234: }
235: }
236: deg = PETSC_INT_MAX;
237: maxDeg = 0;
238: for (PetscInt i = 0; i < Ns; i++) {
239: PetscSpace si;
240: PetscInt iDeg, iMaxDeg;
242: PetscCall(PetscSpaceTensorGetSubspace(sp, i, &si));
243: PetscCall(PetscSpaceGetDegree(si, &iDeg, &iMaxDeg));
244: deg = PetscMin(deg, iDeg);
245: maxDeg += iMaxDeg;
246: }
247: sp->degree = deg;
248: sp->maxDegree = maxDeg;
249: tens->uniform = uniform;
250: tens->setupcalled = PETSC_TRUE;
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: static PetscErrorCode PetscSpaceDestroy_Tensor(PetscSpace sp)
255: {
256: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
257: PetscInt Ns, i;
259: PetscFunctionBegin;
260: Ns = tens->numTensSpaces;
261: if (tens->heightsubspaces) {
262: /* sp->Nv is the spatial dimension, so it is equal to the number
263: * of subspaces on higher co-dimension points */
264: for (PetscInt d = 0; d < sp->Nv; ++d) PetscCall(PetscSpaceDestroy(&tens->heightsubspaces[d]));
265: }
266: PetscCall(PetscFree(tens->heightsubspaces));
267: for (i = 0; i < Ns; i++) PetscCall(PetscSpaceDestroy(&tens->tensspaces[i]));
268: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", NULL));
269: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", NULL));
270: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", NULL));
271: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", NULL));
272: PetscCall(PetscFree(tens->tensspaces));
273: PetscCall(PetscFree(tens));
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: static PetscErrorCode PetscSpaceGetDimension_Tensor(PetscSpace sp, PetscInt *dim)
278: {
279: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
280: PetscInt i, Ns, d;
282: PetscFunctionBegin;
283: PetscCall(PetscSpaceSetUp(sp));
284: Ns = tens->numTensSpaces;
285: d = 1;
286: for (i = 0; i < Ns; i++) {
287: PetscInt id;
289: PetscCall(PetscSpaceGetDimension(tens->tensspaces[i], &id));
290: d *= id;
291: }
292: *dim = d;
293: PetscFunctionReturn(PETSC_SUCCESS);
294: }
296: static PetscErrorCode PetscSpaceEvaluate_Tensor(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
297: {
298: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
299: DM dm = sp->dm;
300: PetscInt Nc = sp->Nc;
301: PetscInt Nv = sp->Nv;
302: PetscInt Ns;
303: PetscReal *lpoints, *sB = NULL, *sD = NULL, *sH = NULL;
304: PetscInt pdim;
306: PetscFunctionBegin;
307: if (!tens->setupcalled) {
308: PetscCall(PetscSpaceSetUp(sp));
309: PetscCall(PetscSpaceEvaluate(sp, npoints, points, B, D, H));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
312: Ns = tens->numTensSpaces;
313: PetscCall(PetscSpaceGetDimension(sp, &pdim));
314: PetscCall(DMGetWorkArray(dm, npoints * Nv, MPIU_REAL, &lpoints));
315: if (B || D || H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &sB));
316: if (D || H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * Nv, MPIU_REAL, &sD));
317: if (H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * Nv * Nv, MPIU_REAL, &sH));
318: if (B) {
319: for (PetscInt i = 0; i < npoints * pdim * Nc; i++) B[i] = 1.;
320: }
321: if (D) {
322: for (PetscInt i = 0; i < npoints * pdim * Nc * Nv; i++) D[i] = 1.;
323: }
324: if (H) {
325: for (PetscInt i = 0; i < npoints * pdim * Nc * Nv * Nv; i++) H[i] = 1.;
326: }
327: for (PetscInt s = 0, d = 0, vstep = 1, cstep = 1; s < Ns; s++) {
328: PetscInt sNv, sNc, spdim;
329: PetscInt vskip, cskip;
331: PetscCall(PetscSpaceGetNumVariables(tens->tensspaces[s], &sNv));
332: PetscCall(PetscSpaceGetNumComponents(tens->tensspaces[s], &sNc));
333: PetscCall(PetscSpaceGetDimension(tens->tensspaces[s], &spdim));
334: PetscCheck(!(pdim % vstep) && !(pdim % spdim), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad tensor loop: Nv %" PetscInt_FMT ", Ns %" PetscInt_FMT ", pdim %" PetscInt_FMT ", s %" PetscInt_FMT ", vstep %" PetscInt_FMT ", spdim %" PetscInt_FMT, Nv, Ns, pdim, s, vstep, spdim);
335: PetscCheck(!(Nc % cstep) && !(Nc % sNc), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad tensor loop: Nv %" PetscInt_FMT ", Ns %" PetscInt_FMT ", Nc %" PetscInt_FMT ", s %" PetscInt_FMT ", cstep %" PetscInt_FMT ", sNc %" PetscInt_FMT, Nv, Ns, Nc, s, cstep, spdim);
336: vskip = pdim / (vstep * spdim);
337: cskip = Nc / (cstep * sNc);
338: for (PetscInt p = 0; p < npoints; ++p) {
339: for (PetscInt i = 0; i < sNv; i++) lpoints[p * sNv + i] = points[p * Nv + d + i];
340: }
341: PetscCall(PetscSpaceEvaluate(tens->tensspaces[s], npoints, lpoints, sB, sD, sH));
342: if (B) {
343: for (PetscInt k = 0; k < vskip; k++) {
344: for (PetscInt si = 0; si < spdim; si++) {
345: for (PetscInt j = 0; j < vstep; j++) {
346: PetscInt i = (k * spdim + si) * vstep + j;
348: for (PetscInt l = 0; l < cskip; l++) {
349: for (PetscInt sc = 0; sc < sNc; sc++) {
350: for (PetscInt m = 0; m < cstep; m++) {
351: PetscInt c = (l * sNc + sc) * cstep + m;
353: for (PetscInt p = 0; p < npoints; p++) B[(pdim * p + i) * Nc + c] *= sB[(spdim * p + si) * sNc + sc];
354: }
355: }
356: }
357: }
358: }
359: }
360: }
361: if (D) {
362: for (PetscInt k = 0; k < vskip; k++) {
363: for (PetscInt si = 0; si < spdim; si++) {
364: for (PetscInt j = 0; j < vstep; j++) {
365: PetscInt i = (k * spdim + si) * vstep + j;
367: for (PetscInt l = 0; l < cskip; l++) {
368: for (PetscInt sc = 0; sc < sNc; sc++) {
369: for (PetscInt m = 0; m < cstep; m++) {
370: PetscInt c = (l * sNc + sc) * cstep + m;
372: for (PetscInt der = 0; der < Nv; der++) {
373: if (der >= d && der < d + sNv) {
374: for (PetscInt p = 0; p < npoints; p++) D[((pdim * p + i) * Nc + c) * Nv + der] *= sD[((spdim * p + si) * sNc + sc) * sNv + der - d];
375: } else {
376: for (PetscInt p = 0; p < npoints; p++) D[((pdim * p + i) * Nc + c) * Nv + der] *= sB[(spdim * p + si) * sNc + sc];
377: }
378: }
379: }
380: }
381: }
382: }
383: }
384: }
385: }
386: if (H) {
387: for (PetscInt k = 0; k < vskip; k++) {
388: for (PetscInt si = 0; si < spdim; si++) {
389: for (PetscInt j = 0; j < vstep; j++) {
390: PetscInt i = (k * spdim + si) * vstep + j;
392: for (PetscInt l = 0; l < cskip; l++) {
393: for (PetscInt sc = 0; sc < sNc; sc++) {
394: for (PetscInt m = 0; m < cstep; m++) {
395: PetscInt c = (l * sNc + sc) * cstep + m;
397: for (PetscInt der = 0; der < Nv; der++) {
398: for (PetscInt der2 = 0; der2 < Nv; der2++) {
399: if (der >= d && der < d + sNv && der2 >= d && der2 < d + sNv) {
400: 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];
401: } else if (der >= d && der < d + sNv) {
402: 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];
403: } else if (der2 >= d && der2 < d + sNv) {
404: 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];
405: } else {
406: for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sB[(spdim * p + si) * sNc + sc];
407: }
408: }
409: }
410: }
411: }
412: }
413: }
414: }
415: }
416: }
417: d += sNv;
418: vstep *= spdim;
419: cstep *= sNc;
420: }
421: if (H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * Nv * Nv, MPIU_REAL, &sH));
422: if (D || H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * Nv, MPIU_REAL, &sD));
423: if (B || D || H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &sB));
424: PetscCall(DMRestoreWorkArray(dm, npoints * Nv, MPIU_REAL, &lpoints));
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: /*@
429: PetscSpaceTensorSetNumSubspaces - Set the number of spaces in the tensor product space
431: Input Parameters:
432: + sp - the function space object
433: - numTensSpaces - the number of spaces
435: Level: intermediate
437: Note:
438: The name NumSubspaces is misleading because it is actually setting the number of defining spaces of the tensor product space, not a number of Subspaces of it
440: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
441: @*/
442: PetscErrorCode PetscSpaceTensorSetNumSubspaces(PetscSpace sp, PetscInt numTensSpaces)
443: {
444: PetscFunctionBegin;
446: PetscTryMethod(sp, "PetscSpaceTensorSetNumSubspaces_C", (PetscSpace, PetscInt), (sp, numTensSpaces));
447: PetscFunctionReturn(PETSC_SUCCESS);
448: }
450: /*@
451: PetscSpaceTensorGetNumSubspaces - Get the number of spaces in the tensor product space
453: Input Parameter:
454: . sp - the function space object
456: Output Parameter:
457: . numTensSpaces - the number of spaces
459: Level: intermediate
461: Note:
462: The name NumSubspaces is misleading because it is actually getting the number of defining spaces of the tensor product space, not a number of Subspaces of it
464: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
465: @*/
466: PetscErrorCode PetscSpaceTensorGetNumSubspaces(PetscSpace sp, PetscInt *numTensSpaces)
467: {
468: PetscFunctionBegin;
470: PetscAssertPointer(numTensSpaces, 2);
471: PetscTryMethod(sp, "PetscSpaceTensorGetNumSubspaces_C", (PetscSpace, PetscInt *), (sp, numTensSpaces));
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: /*@
476: PetscSpaceTensorSetSubspace - Set a space in the tensor product space
478: Input Parameters:
479: + sp - the function space object
480: . s - The space number
481: - subsp - the number of spaces
483: Level: intermediate
485: Note:
486: The name SetSubspace is misleading because it is actually setting one of the defining spaces of the tensor product space, not a Subspace of it
488: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
489: @*/
490: PetscErrorCode PetscSpaceTensorSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
491: {
492: PetscFunctionBegin;
495: PetscTryMethod(sp, "PetscSpaceTensorSetSubspace_C", (PetscSpace, PetscInt, PetscSpace), (sp, s, subsp));
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: /*@
500: PetscSpaceTensorGetSubspace - Get a space in the tensor product space
502: Input Parameters:
503: + sp - the function space object
504: - s - The space number
506: Output Parameter:
507: . subsp - the `PetscSpace`
509: Level: intermediate
511: Note:
512: The name GetSubspace is misleading because it is actually getting one of the defining spaces of the tensor product space, not a Subspace of it
514: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
515: @*/
516: PetscErrorCode PetscSpaceTensorGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
517: {
518: PetscFunctionBegin;
520: PetscAssertPointer(subsp, 3);
521: PetscTryMethod(sp, "PetscSpaceTensorGetSubspace_C", (PetscSpace, PetscInt, PetscSpace *), (sp, s, subsp));
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: static PetscErrorCode PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space, PetscInt numTensSpaces)
526: {
527: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
528: PetscInt Ns;
530: PetscFunctionBegin;
531: PetscCheck(!tens->setupcalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of subspaces after setup called");
532: Ns = tens->numTensSpaces;
533: if (numTensSpaces == Ns) PetscFunctionReturn(PETSC_SUCCESS);
534: if (Ns >= 0) {
535: for (PetscInt s = 0; s < Ns; s++) PetscCall(PetscSpaceDestroy(&tens->tensspaces[s]));
536: PetscCall(PetscFree(tens->tensspaces));
537: }
538: Ns = tens->numTensSpaces = numTensSpaces;
539: PetscCall(PetscCalloc1(Ns, &tens->tensspaces));
540: PetscFunctionReturn(PETSC_SUCCESS);
541: }
543: static PetscErrorCode PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space, PetscInt *numTensSpaces)
544: {
545: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
547: PetscFunctionBegin;
548: *numTensSpaces = tens->numTensSpaces;
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: static PetscErrorCode PetscSpaceTensorSetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace subspace)
553: {
554: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
555: PetscInt Ns;
557: PetscFunctionBegin;
558: PetscCheck(!tens->setupcalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change subspace after setup called");
559: Ns = tens->numTensSpaces;
560: PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceTensorSetNumSubspaces() first");
561: PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
562: PetscCall(PetscObjectReference((PetscObject)subspace));
563: PetscCall(PetscSpaceDestroy(&tens->tensspaces[s]));
564: tens->tensspaces[s] = subspace;
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: static PetscErrorCode PetscSpaceGetHeightSubspace_Tensor(PetscSpace sp, PetscInt height, PetscSpace *subsp)
569: {
570: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
571: PetscInt Nc, dim, order, i;
572: PetscSpace bsp;
574: PetscFunctionBegin;
575: PetscCall(PetscSpaceGetNumVariables(sp, &dim));
576: PetscCheck(tens->uniform && tens->numTensSpaces == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Can only get a generic height subspace of a uniform tensor space of 1d spaces.");
577: PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
578: PetscCall(PetscSpaceGetDegree(sp, &order, NULL));
579: PetscCheck(height <= dim && height >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Asked for space at height %" PetscInt_FMT " for dimension %" PetscInt_FMT " space", height, dim);
580: if (!tens->heightsubspaces) PetscCall(PetscCalloc1(dim, &tens->heightsubspaces));
581: if (height <= dim) {
582: if (!tens->heightsubspaces[height - 1]) {
583: PetscSpace sub;
584: const char *name;
586: PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &bsp));
587: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub));
588: PetscCall(PetscObjectGetName((PetscObject)sp, &name));
589: PetscCall(PetscObjectSetName((PetscObject)sub, name));
590: PetscCall(PetscSpaceSetType(sub, PETSCSPACETENSOR));
591: PetscCall(PetscSpaceSetNumComponents(sub, Nc));
592: PetscCall(PetscSpaceSetDegree(sub, order, PETSC_DETERMINE));
593: PetscCall(PetscSpaceSetNumVariables(sub, dim - height));
594: PetscCall(PetscSpaceTensorSetNumSubspaces(sub, dim - height));
595: for (i = 0; i < dim - height; i++) PetscCall(PetscSpaceTensorSetSubspace(sub, i, bsp));
596: PetscCall(PetscSpaceSetUp(sub));
597: tens->heightsubspaces[height - 1] = sub;
598: }
599: *subsp = tens->heightsubspaces[height - 1];
600: } else {
601: *subsp = NULL;
602: }
603: PetscFunctionReturn(PETSC_SUCCESS);
604: }
606: static PetscErrorCode PetscSpaceTensorGetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace *subspace)
607: {
608: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
609: PetscInt Ns;
611: PetscFunctionBegin;
612: Ns = tens->numTensSpaces;
613: PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceTensorSetNumSubspaces() first");
614: PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
615: *subspace = tens->tensspaces[s];
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: static PetscErrorCode PetscSpaceInitialize_Tensor(PetscSpace sp)
620: {
621: PetscFunctionBegin;
622: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Tensor;
623: sp->ops->setup = PetscSpaceSetUp_Tensor;
624: sp->ops->view = PetscSpaceView_Tensor;
625: sp->ops->destroy = PetscSpaceDestroy_Tensor;
626: sp->ops->getdimension = PetscSpaceGetDimension_Tensor;
627: sp->ops->evaluate = PetscSpaceEvaluate_Tensor;
628: sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Tensor;
629: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", PetscSpaceTensorGetNumSubspaces_Tensor));
630: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", PetscSpaceTensorSetNumSubspaces_Tensor));
631: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", PetscSpaceTensorGetSubspace_Tensor));
632: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", PetscSpaceTensorSetSubspace_Tensor));
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: /*MC
637: PETSCSPACETENSOR = "tensor" - A `PetscSpace` object that encapsulates a tensor product space.
638: A tensor product is created of the components of the subspaces as well.
640: Level: intermediate
642: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`, `PetscSpaceTensorGetSubspace()`, `PetscSpaceTensorSetSubspace()`,
643: `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceTensorSetNumSubspaces()`
644: M*/
646: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Tensor(PetscSpace sp)
647: {
648: PetscSpace_Tensor *tens;
650: PetscFunctionBegin;
652: PetscCall(PetscNew(&tens));
653: sp->data = tens;
655: tens->numTensSpaces = PETSC_DEFAULT;
657: PetscCall(PetscSpaceInitialize_Tensor(sp));
658: PetscFunctionReturn(PETSC_SUCCESS);
659: }