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, "%d_", (int)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 iascii;
124: PetscFunctionBegin;
125: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
126: if (iascii) 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: PetscInt d;
264: /* sp->Nv is the spatial dimension, so it is equal to the number
265: * of subspaces on higher co-dimension points */
266: for (d = 0; d < sp->Nv; ++d) PetscCall(PetscSpaceDestroy(&tens->heightsubspaces[d]));
267: }
268: PetscCall(PetscFree(tens->heightsubspaces));
269: for (i = 0; i < Ns; i++) PetscCall(PetscSpaceDestroy(&tens->tensspaces[i]));
270: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", NULL));
271: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", NULL));
272: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", NULL));
273: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", NULL));
274: PetscCall(PetscFree(tens->tensspaces));
275: PetscCall(PetscFree(tens));
276: PetscFunctionReturn(PETSC_SUCCESS);
277: }
279: static PetscErrorCode PetscSpaceGetDimension_Tensor(PetscSpace sp, PetscInt *dim)
280: {
281: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
282: PetscInt i, Ns, d;
284: PetscFunctionBegin;
285: PetscCall(PetscSpaceSetUp(sp));
286: Ns = tens->numTensSpaces;
287: d = 1;
288: for (i = 0; i < Ns; i++) {
289: PetscInt id;
291: PetscCall(PetscSpaceGetDimension(tens->tensspaces[i], &id));
292: d *= id;
293: }
294: *dim = d;
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: static PetscErrorCode PetscSpaceEvaluate_Tensor(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
299: {
300: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
301: DM dm = sp->dm;
302: PetscInt Nc = sp->Nc;
303: PetscInt Nv = sp->Nv;
304: PetscInt Ns;
305: PetscReal *lpoints, *sB = NULL, *sD = NULL, *sH = NULL;
306: PetscInt pdim;
308: PetscFunctionBegin;
309: if (!tens->setupCalled) {
310: PetscCall(PetscSpaceSetUp(sp));
311: PetscCall(PetscSpaceEvaluate(sp, npoints, points, B, D, H));
312: PetscFunctionReturn(PETSC_SUCCESS);
313: }
314: Ns = tens->numTensSpaces;
315: PetscCall(PetscSpaceGetDimension(sp, &pdim));
316: PetscCall(DMGetWorkArray(dm, npoints * Nv, MPIU_REAL, &lpoints));
317: if (B || D || H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &sB));
318: if (D || H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * Nv, MPIU_REAL, &sD));
319: if (H) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * Nv * Nv, MPIU_REAL, &sH));
320: if (B) {
321: for (PetscInt i = 0; i < npoints * pdim * Nc; i++) B[i] = 1.;
322: }
323: if (D) {
324: for (PetscInt i = 0; i < npoints * pdim * Nc * Nv; i++) D[i] = 1.;
325: }
326: if (H) {
327: for (PetscInt i = 0; i < npoints * pdim * Nc * Nv * Nv; i++) H[i] = 1.;
328: }
329: for (PetscInt s = 0, d = 0, vstep = 1, cstep = 1; s < Ns; s++) {
330: PetscInt sNv, sNc, spdim;
331: PetscInt vskip, cskip;
333: PetscCall(PetscSpaceGetNumVariables(tens->tensspaces[s], &sNv));
334: PetscCall(PetscSpaceGetNumComponents(tens->tensspaces[s], &sNc));
335: PetscCall(PetscSpaceGetDimension(tens->tensspaces[s], &spdim));
336: 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);
337: 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);
338: vskip = pdim / (vstep * spdim);
339: cskip = Nc / (cstep * sNc);
340: for (PetscInt p = 0; p < npoints; ++p) {
341: for (PetscInt i = 0; i < sNv; i++) lpoints[p * sNv + i] = points[p * Nv + d + i];
342: }
343: PetscCall(PetscSpaceEvaluate(tens->tensspaces[s], npoints, lpoints, sB, sD, sH));
344: if (B) {
345: for (PetscInt k = 0; k < vskip; k++) {
346: for (PetscInt si = 0; si < spdim; si++) {
347: for (PetscInt j = 0; j < vstep; j++) {
348: PetscInt i = (k * spdim + si) * vstep + j;
350: for (PetscInt l = 0; l < cskip; l++) {
351: for (PetscInt sc = 0; sc < sNc; sc++) {
352: for (PetscInt m = 0; m < cstep; m++) {
353: PetscInt c = (l * sNc + sc) * cstep + m;
355: for (PetscInt p = 0; p < npoints; p++) B[(pdim * p + i) * Nc + c] *= sB[(spdim * p + si) * sNc + sc];
356: }
357: }
358: }
359: }
360: }
361: }
362: }
363: if (D) {
364: for (PetscInt k = 0; k < vskip; k++) {
365: for (PetscInt si = 0; si < spdim; si++) {
366: for (PetscInt j = 0; j < vstep; j++) {
367: PetscInt i = (k * spdim + si) * vstep + j;
369: for (PetscInt l = 0; l < cskip; l++) {
370: for (PetscInt sc = 0; sc < sNc; sc++) {
371: for (PetscInt m = 0; m < cstep; m++) {
372: PetscInt c = (l * sNc + sc) * cstep + m;
374: for (PetscInt der = 0; der < Nv; der++) {
375: if (der >= d && der < d + sNv) {
376: for (PetscInt p = 0; p < npoints; p++) D[((pdim * p + i) * Nc + c) * Nv + der] *= sD[((spdim * p + si) * sNc + sc) * sNv + der - d];
377: } else {
378: for (PetscInt p = 0; p < npoints; p++) D[((pdim * p + i) * Nc + c) * Nv + der] *= sB[(spdim * p + si) * sNc + sc];
379: }
380: }
381: }
382: }
383: }
384: }
385: }
386: }
387: }
388: if (H) {
389: for (PetscInt k = 0; k < vskip; k++) {
390: for (PetscInt si = 0; si < spdim; si++) {
391: for (PetscInt j = 0; j < vstep; j++) {
392: PetscInt i = (k * spdim + si) * vstep + j;
394: for (PetscInt l = 0; l < cskip; l++) {
395: for (PetscInt sc = 0; sc < sNc; sc++) {
396: for (PetscInt m = 0; m < cstep; m++) {
397: PetscInt c = (l * sNc + sc) * cstep + m;
399: for (PetscInt der = 0; der < Nv; der++) {
400: for (PetscInt der2 = 0; der2 < Nv; der2++) {
401: if (der >= d && der < d + sNv && der2 >= d && der2 < d + sNv) {
402: 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];
403: } else if (der >= d && der < 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 + der - d];
405: } else if (der2 >= d && der2 < d + sNv) {
406: 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];
407: } else {
408: for (PetscInt p = 0; p < npoints; p++) H[(((pdim * p + i) * Nc + c) * Nv + der) * Nv + der2] *= sB[(spdim * p + si) * sNc + sc];
409: }
410: }
411: }
412: }
413: }
414: }
415: }
416: }
417: }
418: }
419: d += sNv;
420: vstep *= spdim;
421: cstep *= sNc;
422: }
423: if (H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * Nv * Nv, MPIU_REAL, &sH));
424: if (D || H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * Nv, MPIU_REAL, &sD));
425: if (B || D || H) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &sB));
426: PetscCall(DMRestoreWorkArray(dm, npoints * Nv, MPIU_REAL, &lpoints));
427: PetscFunctionReturn(PETSC_SUCCESS);
428: }
430: /*@
431: PetscSpaceTensorSetNumSubspaces - Set the number of spaces in the tensor product space
433: Input Parameters:
434: + sp - the function space object
435: - numTensSpaces - the number of spaces
437: Level: intermediate
439: Note:
440: 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
442: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
443: @*/
444: PetscErrorCode PetscSpaceTensorSetNumSubspaces(PetscSpace sp, PetscInt numTensSpaces)
445: {
446: PetscFunctionBegin;
448: PetscTryMethod(sp, "PetscSpaceTensorSetNumSubspaces_C", (PetscSpace, PetscInt), (sp, numTensSpaces));
449: PetscFunctionReturn(PETSC_SUCCESS);
450: }
452: /*@
453: PetscSpaceTensorGetNumSubspaces - Get the number of spaces in the tensor product space
455: Input Parameter:
456: . sp - the function space object
458: Output Parameter:
459: . numTensSpaces - the number of spaces
461: Level: intermediate
463: Note:
464: 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
466: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
467: @*/
468: PetscErrorCode PetscSpaceTensorGetNumSubspaces(PetscSpace sp, PetscInt *numTensSpaces)
469: {
470: PetscFunctionBegin;
472: PetscAssertPointer(numTensSpaces, 2);
473: PetscTryMethod(sp, "PetscSpaceTensorGetNumSubspaces_C", (PetscSpace, PetscInt *), (sp, numTensSpaces));
474: PetscFunctionReturn(PETSC_SUCCESS);
475: }
477: /*@
478: PetscSpaceTensorSetSubspace - Set a space in the tensor product space
480: Input Parameters:
481: + sp - the function space object
482: . s - The space number
483: - subsp - the number of spaces
485: Level: intermediate
487: Note:
488: 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
490: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorGetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
491: @*/
492: PetscErrorCode PetscSpaceTensorSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
493: {
494: PetscFunctionBegin;
497: PetscTryMethod(sp, "PetscSpaceTensorSetSubspace_C", (PetscSpace, PetscInt, PetscSpace), (sp, s, subsp));
498: PetscFunctionReturn(PETSC_SUCCESS);
499: }
501: /*@
502: PetscSpaceTensorGetSubspace - Get a space in the tensor product space
504: Input Parameters:
505: + sp - the function space object
506: - s - The space number
508: Output Parameter:
509: . subsp - the `PetscSpace`
511: Level: intermediate
513: Note:
514: 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
516: .seealso: `PETSCSPACETENSOR`, `PetscSpace`, `PetscSpaceTensorSetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
517: @*/
518: PetscErrorCode PetscSpaceTensorGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
519: {
520: PetscFunctionBegin;
522: PetscAssertPointer(subsp, 3);
523: PetscTryMethod(sp, "PetscSpaceTensorGetSubspace_C", (PetscSpace, PetscInt, PetscSpace *), (sp, s, subsp));
524: PetscFunctionReturn(PETSC_SUCCESS);
525: }
527: static PetscErrorCode PetscSpaceTensorSetNumSubspaces_Tensor(PetscSpace space, PetscInt numTensSpaces)
528: {
529: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
530: PetscInt Ns;
532: PetscFunctionBegin;
533: PetscCheck(!tens->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of subspaces after setup called");
534: Ns = tens->numTensSpaces;
535: if (numTensSpaces == Ns) PetscFunctionReturn(PETSC_SUCCESS);
536: if (Ns >= 0) {
537: PetscInt s;
539: for (s = 0; s < Ns; s++) PetscCall(PetscSpaceDestroy(&tens->tensspaces[s]));
540: PetscCall(PetscFree(tens->tensspaces));
541: }
542: Ns = tens->numTensSpaces = numTensSpaces;
543: PetscCall(PetscCalloc1(Ns, &tens->tensspaces));
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: static PetscErrorCode PetscSpaceTensorGetNumSubspaces_Tensor(PetscSpace space, PetscInt *numTensSpaces)
548: {
549: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
551: PetscFunctionBegin;
552: *numTensSpaces = tens->numTensSpaces;
553: PetscFunctionReturn(PETSC_SUCCESS);
554: }
556: static PetscErrorCode PetscSpaceTensorSetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace subspace)
557: {
558: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
559: PetscInt Ns;
561: PetscFunctionBegin;
562: PetscCheck(!tens->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change subspace after setup called");
563: Ns = tens->numTensSpaces;
564: PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceTensorSetNumSubspaces() first");
565: PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
566: PetscCall(PetscObjectReference((PetscObject)subspace));
567: PetscCall(PetscSpaceDestroy(&tens->tensspaces[s]));
568: tens->tensspaces[s] = subspace;
569: PetscFunctionReturn(PETSC_SUCCESS);
570: }
572: static PetscErrorCode PetscSpaceGetHeightSubspace_Tensor(PetscSpace sp, PetscInt height, PetscSpace *subsp)
573: {
574: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)sp->data;
575: PetscInt Nc, dim, order, i;
576: PetscSpace bsp;
578: PetscFunctionBegin;
579: PetscCall(PetscSpaceGetNumVariables(sp, &dim));
580: 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.");
581: PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
582: PetscCall(PetscSpaceGetDegree(sp, &order, NULL));
583: 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);
584: if (!tens->heightsubspaces) PetscCall(PetscCalloc1(dim, &tens->heightsubspaces));
585: if (height <= dim) {
586: if (!tens->heightsubspaces[height - 1]) {
587: PetscSpace sub;
588: const char *name;
590: PetscCall(PetscSpaceTensorGetSubspace(sp, 0, &bsp));
591: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub));
592: PetscCall(PetscObjectGetName((PetscObject)sp, &name));
593: PetscCall(PetscObjectSetName((PetscObject)sub, name));
594: PetscCall(PetscSpaceSetType(sub, PETSCSPACETENSOR));
595: PetscCall(PetscSpaceSetNumComponents(sub, Nc));
596: PetscCall(PetscSpaceSetDegree(sub, order, PETSC_DETERMINE));
597: PetscCall(PetscSpaceSetNumVariables(sub, dim - height));
598: PetscCall(PetscSpaceTensorSetNumSubspaces(sub, dim - height));
599: for (i = 0; i < dim - height; i++) PetscCall(PetscSpaceTensorSetSubspace(sub, i, bsp));
600: PetscCall(PetscSpaceSetUp(sub));
601: tens->heightsubspaces[height - 1] = sub;
602: }
603: *subsp = tens->heightsubspaces[height - 1];
604: } else {
605: *subsp = NULL;
606: }
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }
610: static PetscErrorCode PetscSpaceTensorGetSubspace_Tensor(PetscSpace space, PetscInt s, PetscSpace *subspace)
611: {
612: PetscSpace_Tensor *tens = (PetscSpace_Tensor *)space->data;
613: PetscInt Ns;
615: PetscFunctionBegin;
616: Ns = tens->numTensSpaces;
617: PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceTensorSetNumSubspaces() first");
618: PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
619: *subspace = tens->tensspaces[s];
620: PetscFunctionReturn(PETSC_SUCCESS);
621: }
623: static PetscErrorCode PetscSpaceInitialize_Tensor(PetscSpace sp)
624: {
625: PetscFunctionBegin;
626: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Tensor;
627: sp->ops->setup = PetscSpaceSetUp_Tensor;
628: sp->ops->view = PetscSpaceView_Tensor;
629: sp->ops->destroy = PetscSpaceDestroy_Tensor;
630: sp->ops->getdimension = PetscSpaceGetDimension_Tensor;
631: sp->ops->evaluate = PetscSpaceEvaluate_Tensor;
632: sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Tensor;
633: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetNumSubspaces_C", PetscSpaceTensorGetNumSubspaces_Tensor));
634: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetNumSubspaces_C", PetscSpaceTensorSetNumSubspaces_Tensor));
635: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorGetSubspace_C", PetscSpaceTensorGetSubspace_Tensor));
636: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceTensorSetSubspace_C", PetscSpaceTensorSetSubspace_Tensor));
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }
640: /*MC
641: PETSCSPACETENSOR = "tensor" - A `PetscSpace` object that encapsulates a tensor product space.
642: A tensor product is created of the components of the subspaces as well.
644: Level: intermediate
646: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`, `PetscSpaceTensorGetSubspace()`, `PetscSpaceTensorSetSubspace()`,
647: `PetscSpaceTensorGetNumSubspaces()`, `PetscSpaceTensorSetNumSubspaces()`
648: M*/
650: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Tensor(PetscSpace sp)
651: {
652: PetscSpace_Tensor *tens;
654: PetscFunctionBegin;
656: PetscCall(PetscNew(&tens));
657: sp->data = tens;
659: tens->numTensSpaces = PETSC_DEFAULT;
661: PetscCall(PetscSpaceInitialize_Tensor(sp));
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }