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;
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.
52: A concatenated sum space will have number of components equal to the sum of the number of components of all subspaces. A non-concatenated,
53: or direct sum space will have the same number of components as its subspaces.
55: Input Parameter:
56: . sp - the function space object
58: Output Parameter:
59: . concatenate - flag indicating whether subspaces are concatenated.
61: Level: intermediate
63: .seealso: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumSetConcatenate()`
64: @*/
65: PetscErrorCode PetscSpaceSumGetConcatenate(PetscSpace sp, PetscBool *concatenate)
66: {
67: PetscFunctionBegin;
69: PetscTryMethod(sp, "PetscSpaceSumGetConcatenate_C", (PetscSpace, PetscBool *), (sp, concatenate));
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: /*@
74: PetscSpaceSumSetConcatenate - Sets the concatenate flag for this space.
75: A concatenated sum space will have number of components equal to the sum of the number of components of all subspaces. A non-concatenated,
76: or direct sum space will have the same number of components as its subspaces .
78: Input Parameters:
79: + sp - the function space object
80: - concatenate - are subspaces concatenated components (true) or direct summands (false)
82: Level: intermediate
84: .seealso: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumGetConcatenate()`
85: @*/
86: PetscErrorCode PetscSpaceSumSetConcatenate(PetscSpace sp, PetscBool concatenate)
87: {
88: PetscFunctionBegin;
90: PetscTryMethod(sp, "PetscSpaceSumSetConcatenate_C", (PetscSpace, PetscBool), (sp, concatenate));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: /*@
95: PetscSpaceSumGetSubspace - Get a space in the sum space
97: Input Parameters:
98: + sp - the function space object
99: - s - The space number
101: Output Parameter:
102: . subsp - the `PetscSpace`
104: Level: intermediate
106: Note:
107: The name GetSubspace is slightly misleading because it is actually getting one of the defining spaces of the sum, not a Subspace of it
109: .seealso: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumSetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
110: @*/
111: PetscErrorCode PetscSpaceSumGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
112: {
113: PetscFunctionBegin;
116: PetscTryMethod(sp, "PetscSpaceSumGetSubspace_C", (PetscSpace, PetscInt, PetscSpace *), (sp, s, subsp));
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: /*@
121: PetscSpaceSumSetSubspace - Set a space in the sum space
123: Input Parameters:
124: + sp - the function space object
125: . s - The space number
126: - subsp - the number of spaces
128: Level: intermediate
130: Note:
131: The name SetSubspace is slightly misleading because it is actually setting one of the defining spaces of the sum, not a Subspace of it
133: .seealso: `PETSCSPACESUM`, `PetscSpace`, `PetscSpaceSumGetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
134: @*/
135: PetscErrorCode PetscSpaceSumSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
136: {
137: PetscFunctionBegin;
140: PetscTryMethod(sp, "PetscSpaceSumSetSubspace_C", (PetscSpace, PetscInt, PetscSpace), (sp, s, subsp));
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: static PetscErrorCode PetscSpaceSumGetNumSubspaces_Sum(PetscSpace space, PetscInt *numSumSpaces)
145: {
146: PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
148: PetscFunctionBegin;
149: *numSumSpaces = sum->numSumSpaces;
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: static PetscErrorCode PetscSpaceSumSetNumSubspaces_Sum(PetscSpace space, PetscInt numSumSpaces)
154: {
155: PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
156: PetscInt Ns = sum->numSumSpaces;
158: PetscFunctionBegin;
159: PetscCheck(!sum->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change number of subspaces after setup called");
160: if (numSumSpaces == Ns) PetscFunctionReturn(PETSC_SUCCESS);
161: if (Ns >= 0) {
162: PetscInt s;
163: for (s = 0; s < Ns; ++s) PetscCall(PetscSpaceDestroy(&sum->sumspaces[s]));
164: PetscCall(PetscFree(sum->sumspaces));
165: }
167: Ns = sum->numSumSpaces = numSumSpaces;
168: PetscCall(PetscCalloc1(Ns, &sum->sumspaces));
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: static PetscErrorCode PetscSpaceSumGetConcatenate_Sum(PetscSpace sp, PetscBool *concatenate)
173: {
174: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
176: PetscFunctionBegin;
177: *concatenate = sum->concatenate;
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: static PetscErrorCode PetscSpaceSumSetConcatenate_Sum(PetscSpace sp, PetscBool concatenate)
182: {
183: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
185: PetscFunctionBegin;
186: PetscCheck(!sum->setupCalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change space concatenation after setup called.");
188: sum->concatenate = concatenate;
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: static PetscErrorCode PetscSpaceSumGetSubspace_Sum(PetscSpace space, PetscInt s, PetscSpace *subspace)
193: {
194: PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
195: PetscInt Ns = sum->numSumSpaces;
197: PetscFunctionBegin;
198: PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceSumSetNumSubspaces() first");
199: PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
201: *subspace = sum->sumspaces[s];
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: static PetscErrorCode PetscSpaceSumSetSubspace_Sum(PetscSpace space, PetscInt s, PetscSpace subspace)
206: {
207: PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
208: PetscInt Ns = sum->numSumSpaces;
210: PetscFunctionBegin;
211: PetscCheck(!sum->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change subspace after setup called");
212: PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscSpaceSumSetNumSubspaces() first");
213: PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);
215: PetscCall(PetscObjectReference((PetscObject)subspace));
216: PetscCall(PetscSpaceDestroy(&sum->sumspaces[s]));
217: sum->sumspaces[s] = subspace;
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: static PetscErrorCode PetscSpaceSetFromOptions_Sum(PetscSpace sp, PetscOptionItems *PetscOptionsObject)
222: {
223: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
224: PetscInt Ns, Nc, Nv, deg, i;
225: PetscBool concatenate = PETSC_TRUE;
226: const char *prefix;
228: PetscFunctionBegin;
229: PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
230: if (!Nv) PetscFunctionReturn(PETSC_SUCCESS);
231: PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
232: PetscCall(PetscSpaceSumGetNumSubspaces(sp, &Ns));
233: PetscCall(PetscSpaceGetDegree(sp, °, NULL));
234: Ns = (Ns == PETSC_DEFAULT) ? 1 : Ns;
236: PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace sum options");
237: PetscCall(PetscOptionsBoundedInt("-petscspace_sum_spaces", "The number of subspaces", "PetscSpaceSumSetNumSubspaces", Ns, &Ns, NULL, 0));
238: PetscCall(PetscOptionsBool("-petscspace_sum_concatenate", "Subspaces are concatenated components of the final space", "PetscSpaceSumSetFromOptions", concatenate, &concatenate, NULL));
239: PetscOptionsHeadEnd();
241: PetscCheck(Ns >= 0 && (Nv <= 0 || Ns != 0), PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have a sum space of %" PetscInt_FMT " spaces", Ns);
242: if (Ns != sum->numSumSpaces) PetscCall(PetscSpaceSumSetNumSubspaces(sp, Ns));
243: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix));
244: for (i = 0; i < Ns; ++i) {
245: PetscInt sNv;
246: PetscSpace subspace;
248: PetscCall(PetscSpaceSumGetSubspace(sp, i, &subspace));
249: if (!subspace) {
250: char subspacePrefix[256];
252: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subspace));
253: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)subspace, prefix));
254: PetscCall(PetscSNPrintf(subspacePrefix, 256, "sumcomp_%" PetscInt_FMT "_", i));
255: PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)subspace, subspacePrefix));
256: } else PetscCall(PetscObjectReference((PetscObject)subspace));
257: PetscCall(PetscSpaceSetFromOptions(subspace));
258: PetscCall(PetscSpaceGetNumVariables(subspace, &sNv));
259: PetscCheck(sNv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Subspace %" PetscInt_FMT " has not been set properly, number of variables is 0.", i);
260: PetscCall(PetscSpaceSumSetSubspace(sp, i, subspace));
261: PetscCall(PetscSpaceDestroy(&subspace));
262: }
263: PetscFunctionReturn(PETSC_SUCCESS);
264: }
266: static PetscErrorCode PetscSpaceSetUp_Sum(PetscSpace sp)
267: {
268: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
269: PetscBool concatenate = PETSC_TRUE;
270: PetscBool uniform;
271: PetscInt Nv, Ns, Nc, i, sum_Nc = 0, deg = PETSC_MAX_INT, maxDeg = PETSC_MIN_INT;
272: PetscInt minNc, maxNc;
274: PetscFunctionBegin;
275: if (sum->setupCalled) PetscFunctionReturn(PETSC_SUCCESS);
277: PetscCall(PetscSpaceGetNumVariables(sp, &Nv));
278: PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
279: PetscCall(PetscSpaceSumGetNumSubspaces(sp, &Ns));
280: if (Ns == PETSC_DEFAULT) {
281: Ns = 1;
282: PetscCall(PetscSpaceSumSetNumSubspaces(sp, Ns));
283: }
284: PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have %" PetscInt_FMT " subspaces", Ns);
285: uniform = PETSC_TRUE;
286: if (Ns) {
287: PetscSpace s0;
289: PetscCall(PetscSpaceSumGetSubspace(sp, 0, &s0));
290: for (PetscInt i = 1; i < Ns; i++) {
291: PetscSpace si;
293: PetscCall(PetscSpaceSumGetSubspace(sp, i, &si));
294: if (si != s0) {
295: uniform = PETSC_FALSE;
296: break;
297: }
298: }
299: }
301: minNc = Nc;
302: maxNc = Nc;
303: for (i = 0; i < Ns; ++i) {
304: PetscInt sNv, sNc, iDeg, iMaxDeg;
305: PetscSpace si;
307: PetscCall(PetscSpaceSumGetSubspace(sp, i, &si));
308: PetscCall(PetscSpaceSetUp(si));
309: PetscCall(PetscSpaceGetNumVariables(si, &sNv));
310: PetscCheck(sNv == Nv, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Subspace %" PetscInt_FMT " has %" PetscInt_FMT " variables, space has %" PetscInt_FMT ".", i, sNv, Nv);
311: PetscCall(PetscSpaceGetNumComponents(si, &sNc));
312: if (i == 0 && sNc == Nc) concatenate = PETSC_FALSE;
313: minNc = PetscMin(minNc, sNc);
314: maxNc = PetscMax(maxNc, sNc);
315: sum_Nc += sNc;
316: PetscCall(PetscSpaceSumGetSubspace(sp, i, &si));
317: PetscCall(PetscSpaceGetDegree(si, &iDeg, &iMaxDeg));
318: deg = PetscMin(deg, iDeg);
319: maxDeg = PetscMax(maxDeg, iMaxDeg);
320: }
322: 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);
323: else PetscCheck(minNc == Nc && maxNc == Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Subspaces must have same number of components as the target space.");
325: sp->degree = deg;
326: sp->maxDegree = maxDeg;
327: sum->concatenate = concatenate;
328: sum->uniform = uniform;
329: sum->setupCalled = PETSC_TRUE;
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: static PetscErrorCode PetscSpaceSumView_Ascii(PetscSpace sp, PetscViewer v)
334: {
335: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
336: PetscBool concatenate = sum->concatenate;
337: PetscInt i, Ns = sum->numSumSpaces;
339: PetscFunctionBegin;
340: if (concatenate) PetscCall(PetscViewerASCIIPrintf(v, "Sum space of %" PetscInt_FMT " concatenated subspaces%s\n", Ns, sum->uniform ? " (all identical)" : ""));
341: else PetscCall(PetscViewerASCIIPrintf(v, "Sum space of %" PetscInt_FMT " subspaces%s\n", Ns, sum->uniform ? " (all identical)" : ""));
342: for (i = 0; i < (sum->uniform ? (Ns > 0 ? 1 : 0) : Ns); ++i) {
343: PetscCall(PetscViewerASCIIPushTab(v));
344: PetscCall(PetscSpaceView(sum->sumspaces[i], v));
345: PetscCall(PetscViewerASCIIPopTab(v));
346: }
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: static PetscErrorCode PetscSpaceView_Sum(PetscSpace sp, PetscViewer viewer)
351: {
352: PetscBool iascii;
354: PetscFunctionBegin;
355: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
356: if (iascii) PetscCall(PetscSpaceSumView_Ascii(sp, viewer));
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: static PetscErrorCode PetscSpaceDestroy_Sum(PetscSpace sp)
361: {
362: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
363: PetscInt i, Ns = sum->numSumSpaces;
365: PetscFunctionBegin;
366: for (i = 0; i < Ns; ++i) PetscCall(PetscSpaceDestroy(&sum->sumspaces[i]));
367: PetscCall(PetscFree(sum->sumspaces));
368: if (sum->heightsubspaces) {
369: PetscInt d;
371: /* sp->Nv is the spatial dimension, so it is equal to the number
372: * of subspaces on higher co-dimension points */
373: for (d = 0; d < sp->Nv; ++d) PetscCall(PetscSpaceDestroy(&sum->heightsubspaces[d]));
374: }
375: PetscCall(PetscFree(sum->heightsubspaces));
376: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetSubspace_C", NULL));
377: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetSubspace_C", NULL));
378: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetNumSubspaces_C", NULL));
379: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetNumSubspaces_C", NULL));
380: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetConcatenate_C", NULL));
381: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetConcatenate_C", NULL));
382: PetscCall(PetscFree(sum));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: static PetscErrorCode PetscSpaceGetDimension_Sum(PetscSpace sp, PetscInt *dim)
387: {
388: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
389: PetscInt i, d = 0, Ns = sum->numSumSpaces;
391: PetscFunctionBegin;
392: if (!sum->setupCalled) {
393: PetscCall(PetscSpaceSetUp(sp));
394: PetscCall(PetscSpaceGetDimension(sp, dim));
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: for (i = 0; i < Ns; ++i) {
399: PetscInt id;
401: PetscCall(PetscSpaceGetDimension(sum->sumspaces[i], &id));
402: d += id;
403: }
405: *dim = d;
406: PetscFunctionReturn(PETSC_SUCCESS);
407: }
409: static PetscErrorCode PetscSpaceEvaluate_Sum(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
410: {
411: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
412: PetscBool concatenate = sum->concatenate;
413: DM dm = sp->dm;
414: PetscInt Nc = sp->Nc, Nv = sp->Nv, Ns = sum->numSumSpaces;
415: PetscInt i, s, offset, ncoffset, pdimfull, numelB, numelD, numelH;
416: PetscReal *sB = NULL, *sD = NULL, *sH = NULL;
418: PetscFunctionBegin;
419: if (!sum->setupCalled) {
420: PetscCall(PetscSpaceSetUp(sp));
421: PetscCall(PetscSpaceEvaluate(sp, npoints, points, B, D, H));
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
424: PetscCall(PetscSpaceGetDimension(sp, &pdimfull));
425: numelB = npoints * pdimfull * Nc;
426: numelD = numelB * Nv;
427: numelH = numelD * Nv;
428: if (B || D || H) PetscCall(DMGetWorkArray(dm, numelB, MPIU_REAL, &sB));
429: if (D || H) PetscCall(DMGetWorkArray(dm, numelD, MPIU_REAL, &sD));
430: if (H) PetscCall(DMGetWorkArray(dm, numelH, MPIU_REAL, &sH));
431: if (B)
432: for (i = 0; i < numelB; ++i) B[i] = 0.;
433: if (D)
434: for (i = 0; i < numelD; ++i) D[i] = 0.;
435: if (H)
436: for (i = 0; i < numelH; ++i) H[i] = 0.;
438: for (s = 0, offset = 0, ncoffset = 0; s < Ns; ++s) {
439: PetscInt sNv, spdim, sNc, p;
441: PetscCall(PetscSpaceGetNumVariables(sum->sumspaces[s], &sNv));
442: PetscCall(PetscSpaceGetNumComponents(sum->sumspaces[s], &sNc));
443: PetscCall(PetscSpaceGetDimension(sum->sumspaces[s], &spdim));
444: PetscCheck(offset + spdim <= pdimfull, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Subspace dimensions exceed target space dimension.");
445: if (s == 0 || !sum->uniform) PetscCall(PetscSpaceEvaluate(sum->sumspaces[s], npoints, points, sB, sD, sH));
446: if (B || D || H) {
447: for (p = 0; p < npoints; ++p) {
448: PetscInt j;
450: for (j = 0; j < spdim; ++j) {
451: PetscInt c;
453: for (c = 0; c < sNc; ++c) {
454: PetscInt compoffset, BInd, sBInd;
456: compoffset = concatenate ? c + ncoffset : c;
457: BInd = (p * pdimfull + j + offset) * Nc + compoffset;
458: sBInd = (p * spdim + j) * sNc + c;
459: if (B) B[BInd] = sB[sBInd];
460: if (D || H) {
461: PetscInt v;
463: for (v = 0; v < Nv; ++v) {
464: PetscInt DInd, sDInd;
466: DInd = BInd * Nv + v;
467: sDInd = sBInd * Nv + v;
468: if (D) D[DInd] = sD[sDInd];
469: if (H) {
470: PetscInt v2;
472: for (v2 = 0; v2 < Nv; ++v2) {
473: PetscInt HInd, sHInd;
475: HInd = DInd * Nv + v2;
476: sHInd = sDInd * Nv + v2;
477: H[HInd] = sH[sHInd];
478: }
479: }
480: }
481: }
482: }
483: }
484: }
485: }
486: offset += spdim;
487: ncoffset += sNc;
488: }
490: if (H) PetscCall(DMRestoreWorkArray(dm, numelH, MPIU_REAL, &sH));
491: if (D || H) PetscCall(DMRestoreWorkArray(dm, numelD, MPIU_REAL, &sD));
492: if (B || D || H) PetscCall(DMRestoreWorkArray(dm, numelB, MPIU_REAL, &sB));
493: PetscFunctionReturn(PETSC_SUCCESS);
494: }
496: static PetscErrorCode PetscSpaceGetHeightSubspace_Sum(PetscSpace sp, PetscInt height, PetscSpace *subsp)
497: {
498: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
499: PetscInt Nc, dim, order;
500: PetscBool tensor;
502: PetscFunctionBegin;
503: PetscCall(PetscSpaceGetNumComponents(sp, &Nc));
504: PetscCall(PetscSpaceGetNumVariables(sp, &dim));
505: PetscCall(PetscSpaceGetDegree(sp, &order, NULL));
506: PetscCall(PetscSpacePolynomialGetTensor(sp, &tensor));
507: 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);
508: if (!sum->heightsubspaces) PetscCall(PetscCalloc1(dim, &sum->heightsubspaces));
509: if (height <= dim) {
510: if (!sum->heightsubspaces[height - 1]) {
511: PetscSpace sub;
512: const char *name;
514: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub));
515: PetscCall(PetscObjectGetName((PetscObject)sp, &name));
516: PetscCall(PetscObjectSetName((PetscObject)sub, name));
517: PetscCall(PetscSpaceSetType(sub, PETSCSPACESUM));
518: PetscCall(PetscSpaceSumSetNumSubspaces(sub, sum->numSumSpaces));
519: PetscCall(PetscSpaceSumSetConcatenate(sub, sum->concatenate));
520: PetscCall(PetscSpaceSetNumComponents(sub, Nc));
521: PetscCall(PetscSpaceSetNumVariables(sub, dim - height));
522: for (PetscInt i = 0; i < sum->numSumSpaces; i++) {
523: PetscSpace subh;
525: PetscCall(PetscSpaceGetHeightSubspace(sum->sumspaces[i], height, &subh));
526: PetscCall(PetscSpaceSumSetSubspace(sub, i, subh));
527: }
528: PetscCall(PetscSpaceSetUp(sub));
529: sum->heightsubspaces[height - 1] = sub;
530: }
531: *subsp = sum->heightsubspaces[height - 1];
532: } else {
533: *subsp = NULL;
534: }
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }
538: static PetscErrorCode PetscSpaceInitialize_Sum(PetscSpace sp)
539: {
540: PetscFunctionBegin;
541: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Sum;
542: sp->ops->setup = PetscSpaceSetUp_Sum;
543: sp->ops->view = PetscSpaceView_Sum;
544: sp->ops->destroy = PetscSpaceDestroy_Sum;
545: sp->ops->getdimension = PetscSpaceGetDimension_Sum;
546: sp->ops->evaluate = PetscSpaceEvaluate_Sum;
547: sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Sum;
549: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetNumSubspaces_C", PetscSpaceSumGetNumSubspaces_Sum));
550: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetNumSubspaces_C", PetscSpaceSumSetNumSubspaces_Sum));
551: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetSubspace_C", PetscSpaceSumGetSubspace_Sum));
552: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetSubspace_C", PetscSpaceSumSetSubspace_Sum));
553: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetConcatenate_C", PetscSpaceSumGetConcatenate_Sum));
554: PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetConcatenate_C", PetscSpaceSumSetConcatenate_Sum));
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }
558: /*MC
559: PETSCSPACESUM = "sum" - A `PetscSpace` object that encapsulates a sum of subspaces.
561: Level: intermediate
563: Note:
564: That sum can either be direct or a concatenation. For example if A and B are spaces each with 2 components,
565: 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
566: same number of variables.
568: .seealso: `PetscSpace`, `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`, `PetscSpaceSumGetNumSubspaces()`, `PetscSpaceSumSetNumSubspaces()`,
569: `PetscSpaceSumGetConcatenate()`, `PetscSpaceSumSetConcatenate()`
570: M*/
571: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Sum(PetscSpace sp)
572: {
573: PetscSpace_Sum *sum;
575: PetscFunctionBegin;
577: PetscCall(PetscNew(&sum));
578: sum->numSumSpaces = PETSC_DEFAULT;
579: sp->data = sum;
580: PetscCall(PetscSpaceInitialize_Sum(sp));
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: PETSC_EXTERN PetscErrorCode PetscSpaceCreateSum(PetscInt numSubspaces, const PetscSpace subspaces[], PetscBool concatenate, PetscSpace *sumSpace)
585: {
586: PetscInt i, Nv, Nc = 0;
588: PetscFunctionBegin;
589: if (sumSpace) PetscCall(PetscSpaceDestroy(sumSpace));
590: PetscCall(PetscSpaceCreate(PetscObjectComm((PetscObject)subspaces[0]), sumSpace));
591: PetscCall(PetscSpaceSetType(*sumSpace, PETSCSPACESUM));
592: PetscCall(PetscSpaceSumSetNumSubspaces(*sumSpace, numSubspaces));
593: PetscCall(PetscSpaceSumSetConcatenate(*sumSpace, concatenate));
594: for (i = 0; i < numSubspaces; ++i) {
595: PetscInt sNc;
597: PetscCall(PetscSpaceSumSetSubspace(*sumSpace, i, subspaces[i]));
598: PetscCall(PetscSpaceGetNumComponents(subspaces[i], &sNc));
599: if (concatenate) Nc += sNc;
600: else Nc = sNc;
601: }
602: PetscCall(PetscSpaceGetNumVariables(subspaces[0], &Nv));
603: PetscCall(PetscSpaceSetNumComponents(*sumSpace, Nc));
604: PetscCall(PetscSpaceSetNumVariables(*sumSpace, Nv));
605: PetscCall(PetscSpaceSetUp(*sumSpace));
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }