Actual source code: spacesum.c
1: #include <petsc/private/petscfeimpl.h>
3: /*@
4: PetscSpaceSumGetNumSubspaces - Get the number of spaces in the sum space
6: Input Parameter:
7: . sp - the function space object
9: Output Parameter:
10: . numSumSpaces - the number of spaces
12: Level: intermediate
14: Note:
15: The name NumSubspaces is slightly misleading because it is actually getting the number of defining spaces of the sum, not a number of Subspaces of it
17: .seealso: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumSetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
18: @*/
19: PetscErrorCode PetscSpaceSumGetNumSubspaces(PetscSpace sp, PetscInt *numSumSpaces)
20: {
21: PetscFunctionBegin;
23: PetscAssertPointer(numSumSpaces, 2);
24: PetscTryMethod(sp, "PetscSpaceSumGetNumSubspaces_C", (PetscSpace, PetscInt *), (sp, numSumSpaces));
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: /*@
29: PetscSpaceSumSetNumSubspaces - Set the number of spaces in the sum space
31: Input Parameters:
32: + sp - the function space object
33: - numSumSpaces - the number of spaces
35: Level: intermediate
37: Note:
38: The name NumSubspaces is slightly misleading because it is actually setting the number of defining spaces of the sum, not a number of Subspaces of it
40: .seealso: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumGetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
41: @*/
42: PetscErrorCode PetscSpaceSumSetNumSubspaces(PetscSpace sp, PetscInt numSumSpaces)
43: {
44: PetscFunctionBegin;
46: PetscTryMethod(sp, "PetscSpaceSumSetNumSubspaces_C", (PetscSpace, PetscInt), (sp, numSumSpaces));
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: /*@
51: PetscSpaceSumGetConcatenate - Get the concatenate flag for this space.
53: Input Parameter:
54: . sp - the function space object
56: Output Parameter:
57: . concatenate - flag indicating whether subspaces are concatenated.
59: Level: intermediate
61: Notes:
62: A concatenated sum space will have the number of components equal to the sum of the number of
63: components of all subspaces. A non-concatenated, or direct sum space will have the same
64: number of components as its subspaces.
66: .seealso: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumSetConcatenate()`
67: @*/
68: PetscErrorCode PetscSpaceSumGetConcatenate(PetscSpace sp, PetscBool *concatenate)
69: {
70: PetscFunctionBegin;
72: PetscTryMethod(sp, "PetscSpaceSumGetConcatenate_C", (PetscSpace, PetscBool *), (sp, concatenate));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: /*@
77: PetscSpaceSumSetConcatenate - Sets the concatenate flag for this space.
79: Input Parameters:
80: + sp - the function space object
81: - concatenate - are subspaces concatenated components (true) or direct summands (false)
83: Level: intermediate
85: Notes:
86: A concatenated sum space will have the number of components equal to the sum of the number of
87: components of all subspaces. A non-concatenated, or direct sum space will have the same
88: number of components as its subspaces .
90: .seealso: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumGetConcatenate()`
91: @*/
92: PetscErrorCode PetscSpaceSumSetConcatenate(PetscSpace sp, PetscBool concatenate)
93: {
94: PetscFunctionBegin;
96: PetscTryMethod(sp, "PetscSpaceSumSetConcatenate_C", (PetscSpace, PetscBool), (sp, concatenate));
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: /*@
101: PetscSpaceSumGetSubspace - Get a space in the sum space
103: Input Parameters:
104: + sp - the function space object
105: - s - The space number
107: Output Parameter:
108: . subsp - the `PetscSpace`
110: Level: intermediate
112: Note:
113: The name GetSubspace is slightly misleading because it is actually getting one of the defining spaces of the sum, not a Subspace of it
115: .seealso: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumSetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
116: @*/
117: PetscErrorCode PetscSpaceSumGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
118: {
119: PetscFunctionBegin;
121: PetscAssertPointer(subsp, 3);
122: PetscTryMethod(sp, "PetscSpaceSumGetSubspace_C", (PetscSpace, PetscInt, PetscSpace *), (sp, s, subsp));
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: /*@
127: PetscSpaceSumSetSubspace - Set a space in the sum space
129: Input Parameters:
130: + sp - the function space object
131: . s - The space number
132: - subsp - the number of spaces
134: Level: intermediate
136: Note:
137: The name SetSubspace is slightly misleading because it is actually setting one of the defining spaces of the sum, not a Subspace of it
139: .seealso: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumGetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
140: @*/
141: PetscErrorCode PetscSpaceSumSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
142: {
143: PetscFunctionBegin;
146: PetscTryMethod(sp, "PetscSpaceSumSetSubspace_C", (PetscSpace, PetscInt, PetscSpace), (sp, s, subsp));
147: PetscFunctionReturn(PETSC_SUCCESS);
148: }
150: static PetscErrorCode PetscSpaceSumGetNumSubspaces_Sum(PetscSpace space, PetscInt *numSumSpaces)
151: {
152: PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
154: PetscFunctionBegin;
155: *numSumSpaces = sum->numSumSpaces;
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: static PetscErrorCode PetscSpaceSumSetNumSubspaces_Sum(PetscSpace space, PetscInt numSumSpaces)
160: {
161: PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
162: PetscInt Ns = sum->numSumSpaces;
164: PetscFunctionBegin;
165: PetscCheck(!sum->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of subspaces after setup called");
166: if (numSumSpaces == Ns) PetscFunctionReturn(PETSC_SUCCESS);
167: if (Ns >= 0) {
168: PetscInt s;
169: for (s = 0; s < Ns; ++s) PetscCall(PetscSpaceDestroy(&sum->sumspaces[s]));
170: PetscCall(PetscFree(sum->sumspaces));
171: }
173: Ns = sum->numSumSpaces = numSumSpaces;
174: PetscCall(PetscCalloc1(Ns, &sum->sumspaces));
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: static PetscErrorCode PetscSpaceSumGetConcatenate_Sum(PetscSpace sp, PetscBool *concatenate)
179: {
180: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
182: PetscFunctionBegin;
183: *concatenate = sum->concatenate;
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: static PetscErrorCode PetscSpaceSumSetConcatenate_Sum(PetscSpace sp, PetscBool concatenate)
188: {
189: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
191: PetscFunctionBegin;
192: PetscCheck(!sum->setupCalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change space concatenation after setup called.");
194: sum->concatenate = concatenate;
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }
198: static PetscErrorCode PetscSpaceSumGetSubspace_Sum(PetscSpace space, PetscInt s, PetscSpace *subspace)
199: {
200: PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
201: PetscInt Ns = sum->numSumSpaces;
203: PetscFunctionBegin;
204: PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceSumSetNumSubspaces() first");
205: PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
207: *subspace = sum->sumspaces[s];
208: PetscFunctionReturn(PETSC_SUCCESS);
209: }
211: static PetscErrorCode PetscSpaceSumSetSubspace_Sum(PetscSpace space, PetscInt s, PetscSpace subspace)
212: {
213: PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
214: PetscInt Ns = sum->numSumSpaces;
216: PetscFunctionBegin;
217: PetscCheck(!sum->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change subspace after setup called");
218: PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceSumSetNumSubspaces() first");
219: PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
221: PetscCall(PetscObjectReference((PetscObject)subspace));
222: PetscCall(PetscSpaceDestroy(&sum->sumspaces[s]));
223: sum->sumspaces[s] = subspace;
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: static PetscErrorCode PetscSpaceSetFromOptions_Sum(PetscSpace sp, PetscOptionItems *PetscOptionsObject)
228: {
229: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
230: PetscInt Ns, Nc, Nv, deg, i;
231: PetscBool concatenate = PETSC_TRUE;
232: const char *prefix;
234: PetscFunctionBegin;
235: PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
236: if (!Nv) PetscFunctionReturn(PETSC_SUCCESS);
237: PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
238: PetscCall(PetscSpaceSumGetNumSubspaces(sp, &Ns));
239: PetscCall(PetscSpaceGetDegree(sp, °, NULL));
240: Ns = (Ns == PETSC_DEFAULT) ? 1 : Ns;
242: PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace sum options");
243: PetscCall(PetscOptionsBoundedInt("-petscspace_sum_spaces", "The number of subspaces", "PetscSpaceSumSetNumSubspaces", Ns, &Ns, NULL, 0));
244: PetscCall(PetscOptionsBool("-petscspace_sum_concatenate", "Subspaces are concatenated components of the final space", "PetscSpaceSumSetFromOptions", concatenate, &concatenate, NULL));
245: PetscOptionsHeadEnd();
247: PetscCheck(Ns >= 0 && (Nv <= 0 || Ns != 0), PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a sum space of %" PetscInt_FMT " spaces", Ns);
248: if (Ns != sum->numSumSpaces) PetscCall(PetscSpaceSumSetNumSubspaces(sp, Ns));
249: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix));
250: for (i = 0; i < Ns; ++i) {
251: PetscInt sNv;
252: PetscSpace subspace;
254: PetscCall(PetscSpaceSumGetSubspace(sp, i, &subspace));
255: if (!subspace) {
256: char subspacePrefix[256];
258: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subspace));
259: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subspace, prefix));
260: PetscCall(PetscSNPrintf(subspacePrefix, 256, "sumcomp_%" PetscInt_FMT "_", i));
261: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subspace, subspacePrefix));
262: } else PetscCall(PetscObjectReference((PetscObject)subspace));
263: PetscCall(PetscSpaceSetFromOptions(subspace));
264: PetscCall(PetscSpaceGetNumVariables(subspace, &sNv));
265: PetscCheck(sNv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Subspace %" PetscInt_FMT " has not been set properly, number of variables is 0.", i);
266: PetscCall(PetscSpaceSumSetSubspace(sp, i, subspace));
267: PetscCall(PetscSpaceDestroy(&subspace));
268: }
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: static PetscErrorCode PetscSpaceSetUp_Sum(PetscSpace sp)
273: {
274: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
275: PetscBool concatenate = PETSC_TRUE;
276: PetscBool uniform;
277: PetscInt Nv, Ns, Nc, i, sum_Nc = 0, deg = PETSC_MAX_INT, maxDeg = PETSC_MIN_INT;
278: PetscInt minNc, maxNc;
280: PetscFunctionBegin;
281: if (sum->setupCalled) PetscFunctionReturn(PETSC_SUCCESS);
283: PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
284: PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
285: PetscCall(PetscSpaceSumGetNumSubspaces(sp, &Ns));
286: if (Ns == PETSC_DEFAULT) {
287: Ns = 1;
288: PetscCall(PetscSpaceSumSetNumSubspaces(sp, Ns));
289: }
290: PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have %" PetscInt_FMT " subspaces", Ns);
291: uniform = PETSC_TRUE;
292: if (Ns) {
293: PetscSpace s0;
295: PetscCall(PetscSpaceSumGetSubspace(sp, 0, &s0));
296: for (PetscInt i = 1; i < Ns; i++) {
297: PetscSpace si;
299: PetscCall(PetscSpaceSumGetSubspace(sp, i, &si));
300: if (si != s0) {
301: uniform = PETSC_FALSE;
302: break;
303: }
304: }
305: }
307: minNc = Nc;
308: maxNc = Nc;
309: for (i = 0; i < Ns; ++i) {
310: PetscInt sNv, sNc, iDeg, iMaxDeg;
311: PetscSpace si;
313: PetscCall(PetscSpaceSumGetSubspace(sp, i, &si));
314: PetscCall(PetscSpaceSetUp(si));
315: PetscCall(PetscSpaceGetNumVariables(si, &sNv));
316: PetscCheck(sNv == Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Subspace %" PetscInt_FMT " has %" PetscInt_FMT " variables, space has %" PetscInt_FMT ".", i, sNv, Nv);
317: PetscCall(PetscSpaceGetNumComponents(si, &sNc));
318: if (i == 0 && sNc == Nc) concatenate = PETSC_FALSE;
319: minNc = PetscMin(minNc, sNc);
320: maxNc = PetscMax(maxNc, sNc);
321: sum_Nc += sNc;
322: PetscCall(PetscSpaceSumGetSubspace(sp, i, &si));
323: PetscCall(PetscSpaceGetDegree(si, &iDeg, &iMaxDeg));
324: deg = PetscMin(deg, iDeg);
325: maxDeg = PetscMax(maxDeg, iMaxDeg);
326: }
328: if (concatenate) PetscCheck(sum_Nc == Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Total number of subspace components (%" PetscInt_FMT ") does not match number of target space components (%" PetscInt_FMT ").", sum_Nc, Nc);
329: else PetscCheck(minNc == Nc && maxNc == Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Subspaces must have same number of components as the target space.");
331: sp->degree = deg;
332: sp->maxDegree = maxDeg;
333: sum->concatenate = concatenate;
334: sum->uniform = uniform;
335: sum->setupCalled = PETSC_TRUE;
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: static PetscErrorCode PetscSpaceSumView_Ascii(PetscSpace sp, PetscViewer v)
340: {
341: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
342: PetscBool concatenate = sum->concatenate;
343: PetscInt i, Ns = sum->numSumSpaces;
345: PetscFunctionBegin;
346: if (concatenate) PetscCall(PetscViewerASCIIPrintf(v, "Sum space of %" PetscInt_FMT " concatenated subspaces%s\n", Ns, sum->uniform ? " (all identical)" : ""));
347: else PetscCall(PetscViewerASCIIPrintf(v, "Sum space of %" PetscInt_FMT " subspaces%s\n", Ns, sum->uniform ? " (all identical)" : ""));
348: for (i = 0; i < (sum->uniform ? (Ns > 0 ? 1 : 0) : Ns); ++i) {
349: PetscCall(PetscViewerASCIIPushTab(v));
350: PetscCall(PetscSpaceView(sum->sumspaces[i], v));
351: PetscCall(PetscViewerASCIIPopTab(v));
352: }
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: static PetscErrorCode PetscSpaceView_Sum(PetscSpace sp, PetscViewer viewer)
357: {
358: PetscBool iascii;
360: PetscFunctionBegin;
361: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
362: if (iascii) PetscCall(PetscSpaceSumView_Ascii(sp, viewer));
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }
366: static PetscErrorCode PetscSpaceDestroy_Sum(PetscSpace sp)
367: {
368: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
369: PetscInt i, Ns = sum->numSumSpaces;
371: PetscFunctionBegin;
372: for (i = 0; i < Ns; ++i) PetscCall(PetscSpaceDestroy(&sum->sumspaces[i]));
373: PetscCall(PetscFree(sum->sumspaces));
374: if (sum->heightsubspaces) {
375: PetscInt d;
377: /* sp->Nv is the spatial dimension, so it is equal to the number
378: * of subspaces on higher co-dimension points */
379: for (d = 0; d < sp->Nv; ++d) PetscCall(PetscSpaceDestroy(&sum->heightsubspaces[d]));
380: }
381: PetscCall(PetscFree(sum->heightsubspaces));
382: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetSubspace_C", NULL));
383: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetSubspace_C", NULL));
384: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetNumSubspaces_C", NULL));
385: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetNumSubspaces_C", NULL));
386: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetConcatenate_C", NULL));
387: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetConcatenate_C", NULL));
388: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetInterleave_C", NULL));
389: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetInterleave_C", NULL));
390: PetscCall(PetscFree(sum));
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: static PetscErrorCode PetscSpaceGetDimension_Sum(PetscSpace sp, PetscInt *dim)
395: {
396: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
397: PetscInt i, d = 0, Ns = sum->numSumSpaces;
399: PetscFunctionBegin;
400: if (!sum->setupCalled) {
401: PetscCall(PetscSpaceSetUp(sp));
402: PetscCall(PetscSpaceGetDimension(sp, dim));
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: for (i = 0; i < Ns; ++i) {
407: PetscInt id;
409: PetscCall(PetscSpaceGetDimension(sum->sumspaces[i], &id));
410: d += id;
411: }
413: *dim = d;
414: PetscFunctionReturn(PETSC_SUCCESS);
415: }
417: static PetscErrorCode PetscSpaceEvaluate_Sum(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
418: {
419: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
420: PetscBool concatenate = sum->concatenate;
421: DM dm = sp->dm;
422: PetscInt Nc = sp->Nc, Nv = sp->Nv, Ns = sum->numSumSpaces;
423: PetscInt i, s, offset, ncoffset, pdimfull, numelB, numelD, numelH;
424: PetscReal *sB = NULL, *sD = NULL, *sH = NULL;
426: PetscFunctionBegin;
427: if (!sum->setupCalled) {
428: PetscCall(PetscSpaceSetUp(sp));
429: PetscCall(PetscSpaceEvaluate(sp, npoints, points, B, D, H));
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
432: PetscCall(PetscSpaceGetDimension(sp, &pdimfull));
433: numelB = npoints * pdimfull * Nc;
434: numelD = numelB * Nv;
435: numelH = numelD * Nv;
436: if (B || D || H) PetscCall(DMGetWorkArray(dm, numelB, MPIU_REAL, &sB));
437: if (D || H) PetscCall(DMGetWorkArray(dm, numelD, MPIU_REAL, &sD));
438: if (H) PetscCall(DMGetWorkArray(dm, numelH, MPIU_REAL, &sH));
439: if (B)
440: for (i = 0; i < numelB; ++i) B[i] = 0.;
441: if (D)
442: for (i = 0; i < numelD; ++i) D[i] = 0.;
443: if (H)
444: for (i = 0; i < numelH; ++i) H[i] = 0.;
446: for (s = 0, offset = 0, ncoffset = 0; s < Ns; ++s) {
447: PetscInt sNv, spdim, sNc, p;
449: PetscCall(PetscSpaceGetNumVariables(sum->sumspaces[s], &sNv));
450: PetscCall(PetscSpaceGetNumComponents(sum->sumspaces[s], &sNc));
451: PetscCall(PetscSpaceGetDimension(sum->sumspaces[s], &spdim));
452: PetscCheck(offset + spdim <= pdimfull, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Subspace dimensions exceed target space dimension.");
453: if (s == 0 || !sum->uniform) PetscCall(PetscSpaceEvaluate(sum->sumspaces[s], npoints, points, sB, sD, sH));
454: if (B || D || H) {
455: for (p = 0; p < npoints; ++p) {
456: PetscInt j;
458: for (j = 0; j < spdim; ++j) {
459: PetscInt c;
460: PetscInt b = sum->interleave_basis ? (j * Ns + s) : (j + offset);
462: for (c = 0; c < sNc; ++c) {
463: PetscInt compoffset, BInd, sBInd;
465: compoffset = concatenate ? (sum->interleave_components ? (c * Ns + s) : (c + ncoffset)) : c;
466: BInd = (p * pdimfull + b) * Nc + compoffset;
467: sBInd = (p * spdim + j) * sNc + c;
468: if (B) B[BInd] = sB[sBInd];
469: if (D || H) {
470: PetscInt v;
472: for (v = 0; v < Nv; ++v) {
473: PetscInt DInd, sDInd;
475: DInd = BInd * Nv + v;
476: sDInd = sBInd * Nv + v;
477: if (D) D[DInd] = sD[sDInd];
478: if (H) {
479: PetscInt v2;
481: for (v2 = 0; v2 < Nv; ++v2) {
482: PetscInt HInd, sHInd;
484: HInd = DInd * Nv + v2;
485: sHInd = sDInd * Nv + v2;
486: H[HInd] = sH[sHInd];
487: }
488: }
489: }
490: }
491: }
492: }
493: }
494: }
495: offset += spdim;
496: ncoffset += sNc;
497: }
499: if (H) PetscCall(DMRestoreWorkArray(dm, numelH, MPIU_REAL, &sH));
500: if (D || H) PetscCall(DMRestoreWorkArray(dm, numelD, MPIU_REAL, &sD));
501: if (B || D || H) PetscCall(DMRestoreWorkArray(dm, numelB, MPIU_REAL, &sB));
502: PetscFunctionReturn(PETSC_SUCCESS);
503: }
505: static PetscErrorCode PetscSpaceGetHeightSubspace_Sum(PetscSpace sp, PetscInt height, PetscSpace *subsp)
506: {
507: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
508: PetscInt Nc, dim, order;
509: PetscBool tensor;
511: PetscFunctionBegin;
512: PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
513: PetscCall(PetscSpaceGetNumVariables(sp, &dim));
514: PetscCall(PetscSpaceGetDegree(sp, &order, NULL));
515: PetscCall(PetscSpacePolynomialGetTensor(sp, &tensor));
516: 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);
517: if (!sum->heightsubspaces) PetscCall(PetscCalloc1(dim, &sum->heightsubspaces));
518: if (height <= dim) {
519: if (!sum->heightsubspaces[height - 1]) {
520: PetscSpace sub;
521: const char *name;
523: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub));
524: PetscCall(PetscObjectGetName((PetscObject)sp, &name));
525: PetscCall(PetscObjectSetName((PetscObject)sub, name));
526: PetscCall(PetscSpaceSetType(sub, PETSCSPACESUM));
527: PetscCall(PetscSpaceSumSetNumSubspaces(sub, sum->numSumSpaces));
528: PetscCall(PetscSpaceSumSetConcatenate(sub, sum->concatenate));
529: PetscCall(PetscSpaceSetNumComponents(sub, Nc));
530: PetscCall(PetscSpaceSetNumVariables(sub, dim - height));
531: for (PetscInt i = 0; i < sum->numSumSpaces; i++) {
532: PetscSpace subh;
534: PetscCall(PetscSpaceGetHeightSubspace(sum->sumspaces[i], height, &subh));
535: PetscCall(PetscSpaceSumSetSubspace(sub, i, subh));
536: }
537: PetscCall(PetscSpaceSetUp(sub));
538: sum->heightsubspaces[height - 1] = sub;
539: }
540: *subsp = sum->heightsubspaces[height - 1];
541: } else {
542: *subsp = NULL;
543: }
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: /*@
548: PetscSpaceSumSetInterleave - Set whether the basis functions and components of a uniform sum are interleaved
550: Logically collective
552: Input Parameters:
553: + sp - a `PetscSpace` of type `PETSCSPACESUM`
554: . interleave_basis - if `PETSC_TRUE`, the basis vectors of the subspaces are interleaved
555: - interleave_components - if `PETSC_TRUE` and the space concatenates components (`PetscSpaceSumGetConcatenate()`),
556: interleave the concatenated components
558: Level: developer
560: .seealso: `PetscSpace`, `PETSCSPACESUM`, `PETSCFEVECTOR`, `PetscSpaceSumGetInterleave()`
561: @*/
562: PetscErrorCode PetscSpaceSumSetInterleave(PetscSpace sp, PetscBool interleave_basis, PetscBool interleave_components)
563: {
564: PetscFunctionBegin;
566: PetscTryMethod(sp, "PetscSpaceSumSetInterleave_C", (PetscSpace, PetscBool, PetscBool), (sp, interleave_basis, interleave_components));
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: static PetscErrorCode PetscSpaceSumSetInterleave_Sum(PetscSpace sp, PetscBool interleave_basis, PetscBool interleave_components)
571: {
572: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
574: PetscFunctionBegin;
575: sum->interleave_basis = interleave_basis;
576: sum->interleave_components = interleave_components;
577: PetscFunctionReturn(PETSC_SUCCESS);
578: }
580: /*@
581: PetscSpaceSumGetInterleave - Get whether the basis functions and components of a uniform sum are interleaved
583: Logically collective
585: Input Parameter:
586: . sp - a `PetscSpace` of type `PETSCSPACESUM`
588: Output Parameters:
589: + interleave_basis - if `PETSC_TRUE`, the basis vectors of the subspaces are interleaved
590: - interleave_components - if `PETSC_TRUE` and the space concatenates components (`PetscSpaceSumGetConcatenate()`),
591: interleave the concatenated components
593: Level: developer
595: .seealso: `PetscSpace`, `PETSCSPACESUM`, `PETSCFEVECTOR`, `PetscSpaceSumSetInterleave()`
596: @*/
597: PetscErrorCode PetscSpaceSumGetInterleave(PetscSpace sp, PetscBool *interleave_basis, PetscBool *interleave_components)
598: {
599: PetscFunctionBegin;
601: if (interleave_basis) PetscAssertPointer(interleave_basis, 2);
602: if (interleave_components) PetscAssertPointer(interleave_components, 3);
603: PetscTryMethod(sp, "PetscSpaceSumGetInterleave_C", (PetscSpace, PetscBool *, PetscBool *), (sp, interleave_basis, interleave_components));
604: PetscFunctionReturn(PETSC_SUCCESS);
605: }
607: static PetscErrorCode PetscSpaceSumGetInterleave_Sum(PetscSpace sp, PetscBool *interleave_basis, PetscBool *interleave_components)
608: {
609: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
611: PetscFunctionBegin;
612: if (interleave_basis) *interleave_basis = sum->interleave_basis;
613: if (interleave_components) *interleave_components = sum->interleave_components;
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: static PetscErrorCode PetscSpaceInitialize_Sum(PetscSpace sp)
618: {
619: PetscFunctionBegin;
620: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Sum;
621: sp->ops->setup = PetscSpaceSetUp_Sum;
622: sp->ops->view = PetscSpaceView_Sum;
623: sp->ops->destroy = PetscSpaceDestroy_Sum;
624: sp->ops->getdimension = PetscSpaceGetDimension_Sum;
625: sp->ops->evaluate = PetscSpaceEvaluate_Sum;
626: sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Sum;
628: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetNumSubspaces_C", PetscSpaceSumGetNumSubspaces_Sum));
629: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetNumSubspaces_C", PetscSpaceSumSetNumSubspaces_Sum));
630: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetSubspace_C", PetscSpaceSumGetSubspace_Sum));
631: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetSubspace_C", PetscSpaceSumSetSubspace_Sum));
632: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetConcatenate_C", PetscSpaceSumGetConcatenate_Sum));
633: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetConcatenate_C", PetscSpaceSumSetConcatenate_Sum));
634: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetInterleave_C", PetscSpaceSumGetInterleave_Sum));
635: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetInterleave_C", PetscSpaceSumSetInterleave_Sum));
636: PetscFunctionReturn(PETSC_SUCCESS);
637: }
639: /*MC
640: PETSCSPACESUM = "sum" - A `PetscSpace` object that encapsulates a sum of subspaces.
642: Level: intermediate
644: Note:
645: That sum can either be direct or a concatenation. For example if A and B are spaces each with 2 components,
646: the direct sum of A and B will also have 2 components while the concatenated sum will have 4 components. In both cases A and B must be defined over the
647: same number of variables.
649: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`, `PetscSpaceSumGetNumSubspaces()`, `PetscSpaceSumSetNumSubspaces()`,
650: `PetscSpaceSumGetConcatenate()`, `PetscSpaceSumSetConcatenate()`
651: M*/
652: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Sum(PetscSpace sp)
653: {
654: PetscSpace_Sum *sum;
656: PetscFunctionBegin;
658: PetscCall(PetscNew(&sum));
659: sum->numSumSpaces = PETSC_DEFAULT;
660: sp->data = sum;
661: PetscCall(PetscSpaceInitialize_Sum(sp));
662: PetscFunctionReturn(PETSC_SUCCESS);
663: }
665: PETSC_EXTERN PetscErrorCode PetscSpaceCreateSum(PetscInt numSubspaces, const PetscSpace subspaces[], PetscBool concatenate, PetscSpace *sumSpace)
666: {
667: PetscInt i, Nv, Nc = 0;
669: PetscFunctionBegin;
670: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)subspaces[0]), sumSpace));
671: PetscCall(PetscSpaceSetType(*sumSpace, PETSCSPACESUM));
672: PetscCall(PetscSpaceSumSetNumSubspaces(*sumSpace, numSubspaces));
673: PetscCall(PetscSpaceSumSetConcatenate(*sumSpace, concatenate));
674: for (i = 0; i < numSubspaces; ++i) {
675: PetscInt sNc;
677: PetscCall(PetscSpaceSumSetSubspace(*sumSpace, i, subspaces[i]));
678: PetscCall(PetscSpaceGetNumComponents(subspaces[i], &sNc));
679: if (concatenate) Nc += sNc;
680: else Nc = sNc;
681: }
682: PetscCall(PetscSpaceGetNumVariables(subspaces[0], &Nv));
683: PetscCall(PetscSpaceSetNumComponents(*sumSpace, Nc));
684: PetscCall(PetscSpaceSetNumVariables(*sumSpace, Nv));
685: PetscCall(PetscSpaceSetUp(*sumSpace));
686: PetscFunctionReturn(PETSC_SUCCESS);
687: }