Actual source code: spacesum.c
1: #include <petsc/private/petscfeimpl.h>
2: /*@
3: PetscSpaceSumGetNumSubspaces - Get the number of spaces in the sum
5: Input Parameter:
6: . sp - the function space object
8: Output Parameter:
9: . numSumSpaces - the number of spaces
11: Level: intermediate
13: .seealso: `PetscSpaceSumSetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
14: @*/
15: PetscErrorCode PetscSpaceSumGetNumSubspaces(PetscSpace sp, PetscInt *numSumSpaces)
16: {
19: PetscTryMethod(sp, "PetscSpaceSumGetNumSubspaces_C", (PetscSpace, PetscInt *), (sp, numSumSpaces));
20: return 0;
21: }
23: /*@
24: PetscSpaceSumSetNumSubspaces - Set the number of spaces in the sum
26: Input Parameters:
27: + sp - the function space object
28: - numSumSpaces - the number of spaces
30: Level: intermediate
32: .seealso: `PetscSpaceSumGetNumSubspaces()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
33: @*/
34: PetscErrorCode PetscSpaceSumSetNumSubspaces(PetscSpace sp, PetscInt numSumSpaces)
35: {
37: PetscTryMethod(sp, "PetscSpaceSumSetNumSubspaces_C", (PetscSpace, PetscInt), (sp, numSumSpaces));
38: return 0;
39: }
41: /*@
42: PetscSpaceSumGetConcatenate - Get the concatenate flag for this space.
43: A concatenated sum space will have number of components equal to the sum of the number of components of all subspaces.A non-concatenated,
44: or direct sum space will have the same number of components as its subspaces .
46: Input Parameters:
47: . sp - the function space object
49: Output Parameters:
50: . concatenate - flag indicating whether subspaces are concatenated.
52: Level: intermediate
54: .seealso: `PetscSpaceSumSetConcatenate()`
55: @*/
56: PetscErrorCode PetscSpaceSumGetConcatenate(PetscSpace sp, PetscBool *concatenate)
57: {
59: PetscTryMethod(sp, "PetscSpaceSumGetConcatenate_C", (PetscSpace, PetscBool *), (sp, concatenate));
60: return 0;
61: }
63: /*@
64: PetscSpaceSumSetConcatenate - Sets the concatenate flag for this space.
65: A concatenated sum space will have number of components equal to the sum of the number of components of all subspaces.A non-concatenated,
66: or direct sum space will have the same number of components as its subspaces .
68: Input Parameters:
69: + sp - the function space object
70: - concatenate - are subspaces concatenated components (true) or direct summands (false)
72: Level: intermediate
73: .seealso: `PetscSpaceSumGetConcatenate()`
74: @*/
75: PetscErrorCode PetscSpaceSumSetConcatenate(PetscSpace sp, PetscBool concatenate)
76: {
78: PetscTryMethod(sp, "PetscSpaceSumSetConcatenate_C", (PetscSpace, PetscBool), (sp, concatenate));
79: return 0;
80: }
82: /*@
83: PetscSpaceSumGetSubspace - Get a space in the sum
85: Input Parameters:
86: + sp - the function space object
87: - s - The space number
89: Output Parameter:
90: . subsp - the PetscSpace
92: Level: intermediate
94: .seealso: `PetscSpaceSumSetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
95: @*/
96: PetscErrorCode PetscSpaceSumGetSubspace(PetscSpace sp, PetscInt s, PetscSpace *subsp)
97: {
100: PetscTryMethod(sp, "PetscSpaceSumGetSubspace_C", (PetscSpace, PetscInt, PetscSpace *), (sp, s, subsp));
101: return 0;
102: }
104: /*@
105: PetscSpaceSumSetSubspace - Set a space in the sum
107: Input Parameters:
108: + sp - the function space object
109: . s - The space number
110: - subsp - the number of spaces
112: Level: intermediate
114: .seealso: `PetscSpaceSumGetSubspace()`, `PetscSpaceSetDegree()`, `PetscSpaceSetNumVariables()`
115: @*/
116: PetscErrorCode PetscSpaceSumSetSubspace(PetscSpace sp, PetscInt s, PetscSpace subsp)
117: {
120: PetscTryMethod(sp, "PetscSpaceSumSetSubspace_C", (PetscSpace, PetscInt, PetscSpace), (sp, s, subsp));
121: return 0;
122: }
124: static PetscErrorCode PetscSpaceSumGetNumSubspaces_Sum(PetscSpace space, PetscInt *numSumSpaces)
125: {
126: PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
128: *numSumSpaces = sum->numSumSpaces;
129: return 0;
130: }
132: static PetscErrorCode PetscSpaceSumSetNumSubspaces_Sum(PetscSpace space, PetscInt numSumSpaces)
133: {
134: PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
135: PetscInt Ns = sum->numSumSpaces;
138: if (numSumSpaces == Ns) return 0;
139: if (Ns >= 0) {
140: PetscInt s;
141: for (s = 0; s < Ns; ++s) PetscSpaceDestroy(&sum->sumspaces[s]);
142: PetscFree(sum->sumspaces);
143: }
145: Ns = sum->numSumSpaces = numSumSpaces;
146: PetscCalloc1(Ns, &sum->sumspaces);
147: return 0;
148: }
150: static PetscErrorCode PetscSpaceSumGetConcatenate_Sum(PetscSpace sp, PetscBool *concatenate)
151: {
152: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
154: *concatenate = sum->concatenate;
155: return 0;
156: }
158: static PetscErrorCode PetscSpaceSumSetConcatenate_Sum(PetscSpace sp, PetscBool concatenate)
159: {
160: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
164: sum->concatenate = concatenate;
165: return 0;
166: }
168: static PetscErrorCode PetscSpaceSumGetSubspace_Sum(PetscSpace space, PetscInt s, PetscSpace *subspace)
169: {
170: PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
171: PetscInt Ns = sum->numSumSpaces;
176: *subspace = sum->sumspaces[s];
177: return 0;
178: }
180: static PetscErrorCode PetscSpaceSumSetSubspace_Sum(PetscSpace space, PetscInt s, PetscSpace subspace)
181: {
182: PetscSpace_Sum *sum = (PetscSpace_Sum *)space->data;
183: PetscInt Ns = sum->numSumSpaces;
189: PetscObjectReference((PetscObject)subspace);
190: PetscSpaceDestroy(&sum->sumspaces[s]);
191: sum->sumspaces[s] = subspace;
192: return 0;
193: }
195: static PetscErrorCode PetscSpaceSetFromOptions_Sum(PetscSpace sp, PetscOptionItems *PetscOptionsObject)
196: {
197: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
198: PetscInt Ns, Nc, Nv, deg, i;
199: PetscBool concatenate = PETSC_TRUE;
200: const char *prefix;
202: PetscSpaceGetNumVariables(sp, &Nv);
203: if (!Nv) return 0;
204: PetscSpaceGetNumComponents(sp, &Nc);
205: PetscSpaceSumGetNumSubspaces(sp, &Ns);
206: PetscSpaceGetDegree(sp, °, NULL);
207: Ns = (Ns == PETSC_DEFAULT) ? 1 : Ns;
209: PetscOptionsHeadBegin(PetscOptionsObject, "PetscSpace sum options");
210: PetscOptionsBoundedInt("-petscspace_sum_spaces", "The number of subspaces", "PetscSpaceSumSetNumSubspaces", Ns, &Ns, NULL, 0);
211: PetscOptionsBool("-petscspace_sum_concatenate", "Subspaces are concatenated components of the final space", "PetscSpaceSumSetFromOptions", concatenate, &concatenate, NULL);
212: PetscOptionsHeadEnd();
215: if (Ns != sum->numSumSpaces) PetscSpaceSumSetNumSubspaces(sp, Ns);
216: PetscObjectGetOptionsPrefix((PetscObject)sp, &prefix);
217: for (i = 0; i < Ns; ++i) {
218: PetscInt sNv;
219: PetscSpace subspace;
221: PetscSpaceSumGetSubspace(sp, i, &subspace);
222: if (!subspace) {
223: char subspacePrefix[256];
225: PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &subspace);
226: PetscObjectSetOptionsPrefix((PetscObject)subspace, prefix);
227: PetscSNPrintf(subspacePrefix, 256, "sumcomp_%" PetscInt_FMT "_", i);
228: PetscObjectAppendOptionsPrefix((PetscObject)subspace, subspacePrefix);
229: } else PetscObjectReference((PetscObject)subspace);
230: PetscSpaceSetFromOptions(subspace);
231: PetscSpaceGetNumVariables(subspace, &sNv);
233: PetscSpaceSumSetSubspace(sp, i, subspace);
234: PetscSpaceDestroy(&subspace);
235: }
236: return 0;
237: }
239: static PetscErrorCode PetscSpaceSetUp_Sum(PetscSpace sp)
240: {
241: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
242: PetscBool concatenate = PETSC_TRUE;
243: PetscBool uniform;
244: PetscInt Nv, Ns, Nc, i, sum_Nc = 0, deg = PETSC_MAX_INT, maxDeg = PETSC_MIN_INT;
245: PetscInt minNc, maxNc;
247: if (sum->setupCalled) return 0;
249: PetscSpaceGetNumVariables(sp, &Nv);
250: PetscSpaceGetNumComponents(sp, &Nc);
251: PetscSpaceSumGetNumSubspaces(sp, &Ns);
252: if (Ns == PETSC_DEFAULT) {
253: Ns = 1;
254: PetscSpaceSumSetNumSubspaces(sp, Ns);
255: }
257: uniform = PETSC_TRUE;
258: if (Ns) {
259: PetscSpace s0;
261: PetscSpaceSumGetSubspace(sp, 0, &s0);
262: for (PetscInt i = 1; i < Ns; i++) {
263: PetscSpace si;
265: PetscSpaceSumGetSubspace(sp, i, &si);
266: if (si != s0) {
267: uniform = PETSC_FALSE;
268: break;
269: }
270: }
271: }
273: minNc = Nc;
274: maxNc = Nc;
275: for (i = 0; i < Ns; ++i) {
276: PetscInt sNv, sNc, iDeg, iMaxDeg;
277: PetscSpace si;
279: PetscSpaceSumGetSubspace(sp, i, &si);
280: PetscSpaceSetUp(si);
281: PetscSpaceGetNumVariables(si, &sNv);
283: PetscSpaceGetNumComponents(si, &sNc);
284: if (i == 0 && sNc == Nc) concatenate = PETSC_FALSE;
285: minNc = PetscMin(minNc, sNc);
286: maxNc = PetscMax(maxNc, sNc);
287: sum_Nc += sNc;
288: PetscSpaceSumGetSubspace(sp, i, &si);
289: PetscSpaceGetDegree(si, &iDeg, &iMaxDeg);
290: deg = PetscMin(deg, iDeg);
291: maxDeg = PetscMax(maxDeg, iMaxDeg);
292: }
297: sp->degree = deg;
298: sp->maxDegree = maxDeg;
299: sum->concatenate = concatenate;
300: sum->uniform = uniform;
301: sum->setupCalled = PETSC_TRUE;
302: return 0;
303: }
305: static PetscErrorCode PetscSpaceSumView_Ascii(PetscSpace sp, PetscViewer v)
306: {
307: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
308: PetscBool concatenate = sum->concatenate;
309: PetscInt i, Ns = sum->numSumSpaces;
311: if (concatenate) PetscViewerASCIIPrintf(v, "Sum space of %" PetscInt_FMT " concatenated subspaces%s\n", Ns, sum->uniform ? " (all identical)" : "");
312: else PetscViewerASCIIPrintf(v, "Sum space of %" PetscInt_FMT " subspaces%s\n", Ns, sum->uniform ? " (all identical)" : "");
313: for (i = 0; i < (sum->uniform ? (Ns > 0 ? 1 : 0) : Ns); ++i) {
314: PetscViewerASCIIPushTab(v);
315: PetscSpaceView(sum->sumspaces[i], v);
316: PetscViewerASCIIPopTab(v);
317: }
318: return 0;
319: }
321: static PetscErrorCode PetscSpaceView_Sum(PetscSpace sp, PetscViewer viewer)
322: {
323: PetscBool iascii;
325: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
326: if (iascii) PetscSpaceSumView_Ascii(sp, viewer);
327: return 0;
328: }
330: static PetscErrorCode PetscSpaceDestroy_Sum(PetscSpace sp)
331: {
332: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
333: PetscInt i, Ns = sum->numSumSpaces;
335: for (i = 0; i < Ns; ++i) PetscSpaceDestroy(&sum->sumspaces[i]);
336: PetscFree(sum->sumspaces);
337: if (sum->heightsubspaces) {
338: PetscInt d;
340: /* sp->Nv is the spatial dimension, so it is equal to the number
341: * of subspaces on higher co-dimension points */
342: for (d = 0; d < sp->Nv; ++d) PetscSpaceDestroy(&sum->heightsubspaces[d]);
343: }
344: PetscFree(sum->heightsubspaces);
345: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetSubspace_C", NULL);
346: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetSubspace_C", NULL);
347: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetNumSubspaces_C", NULL);
348: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetNumSubspaces_C", NULL);
349: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetConcatenate_C", NULL);
350: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetConcatenate_C", NULL);
351: PetscFree(sum);
352: return 0;
353: }
355: static PetscErrorCode PetscSpaceGetDimension_Sum(PetscSpace sp, PetscInt *dim)
356: {
357: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
358: PetscInt i, d = 0, Ns = sum->numSumSpaces;
360: if (!sum->setupCalled) {
361: PetscSpaceSetUp(sp);
362: PetscSpaceGetDimension(sp, dim);
363: return 0;
364: }
366: for (i = 0; i < Ns; ++i) {
367: PetscInt id;
369: PetscSpaceGetDimension(sum->sumspaces[i], &id);
370: d += id;
371: }
373: *dim = d;
374: return 0;
375: }
377: static PetscErrorCode PetscSpaceEvaluate_Sum(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
378: {
379: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
380: PetscBool concatenate = sum->concatenate;
381: DM dm = sp->dm;
382: PetscInt Nc = sp->Nc, Nv = sp->Nv, Ns = sum->numSumSpaces;
383: PetscInt i, s, offset, ncoffset, pdimfull, numelB, numelD, numelH;
384: PetscReal *sB = NULL, *sD = NULL, *sH = NULL;
386: if (!sum->setupCalled) {
387: PetscSpaceSetUp(sp);
388: PetscSpaceEvaluate(sp, npoints, points, B, D, H);
389: return 0;
390: }
391: PetscSpaceGetDimension(sp, &pdimfull);
392: numelB = npoints * pdimfull * Nc;
393: numelD = numelB * Nv;
394: numelH = numelD * Nv;
395: if (B || D || H) DMGetWorkArray(dm, numelB, MPIU_REAL, &sB);
396: if (D || H) DMGetWorkArray(dm, numelD, MPIU_REAL, &sD);
397: if (H) DMGetWorkArray(dm, numelH, MPIU_REAL, &sH);
398: if (B)
399: for (i = 0; i < numelB; ++i) B[i] = 0.;
400: if (D)
401: for (i = 0; i < numelD; ++i) D[i] = 0.;
402: if (H)
403: for (i = 0; i < numelH; ++i) H[i] = 0.;
405: for (s = 0, offset = 0, ncoffset = 0; s < Ns; ++s) {
406: PetscInt sNv, spdim, sNc, p;
408: PetscSpaceGetNumVariables(sum->sumspaces[s], &sNv);
409: PetscSpaceGetNumComponents(sum->sumspaces[s], &sNc);
410: PetscSpaceGetDimension(sum->sumspaces[s], &spdim);
412: if (s == 0 || !sum->uniform) PetscSpaceEvaluate(sum->sumspaces[s], npoints, points, sB, sD, sH);
413: if (B || D || H) {
414: for (p = 0; p < npoints; ++p) {
415: PetscInt j;
417: for (j = 0; j < spdim; ++j) {
418: PetscInt c;
420: for (c = 0; c < sNc; ++c) {
421: PetscInt compoffset, BInd, sBInd;
423: compoffset = concatenate ? c + ncoffset : c;
424: BInd = (p * pdimfull + j + offset) * Nc + compoffset;
425: sBInd = (p * spdim + j) * sNc + c;
426: if (B) B[BInd] = sB[sBInd];
427: if (D || H) {
428: PetscInt v;
430: for (v = 0; v < Nv; ++v) {
431: PetscInt DInd, sDInd;
433: DInd = BInd * Nv + v;
434: sDInd = sBInd * Nv + v;
435: if (D) D[DInd] = sD[sDInd];
436: if (H) {
437: PetscInt v2;
439: for (v2 = 0; v2 < Nv; ++v2) {
440: PetscInt HInd, sHInd;
442: HInd = DInd * Nv + v2;
443: sHInd = sDInd * Nv + v2;
444: H[HInd] = sH[sHInd];
445: }
446: }
447: }
448: }
449: }
450: }
451: }
452: }
453: offset += spdim;
454: ncoffset += sNc;
455: }
457: if (H) DMRestoreWorkArray(dm, numelH, MPIU_REAL, &sH);
458: if (D || H) DMRestoreWorkArray(dm, numelD, MPIU_REAL, &sD);
459: if (B || D || H) DMRestoreWorkArray(dm, numelB, MPIU_REAL, &sB);
460: return 0;
461: }
463: static PetscErrorCode PetscSpaceGetHeightSubspace_Sum(PetscSpace sp, PetscInt height, PetscSpace *subsp)
464: {
465: PetscSpace_Sum *sum = (PetscSpace_Sum *)sp->data;
466: PetscInt Nc, dim, order;
467: PetscBool tensor;
469: PetscSpaceGetNumComponents(sp, &Nc);
470: PetscSpaceGetNumVariables(sp, &dim);
471: PetscSpaceGetDegree(sp, &order, NULL);
472: PetscSpacePolynomialGetTensor(sp, &tensor);
474: if (!sum->heightsubspaces) PetscCalloc1(dim, &sum->heightsubspaces);
475: if (height <= dim) {
476: if (!sum->heightsubspaces[height - 1]) {
477: PetscSpace sub;
478: const char *name;
480: PetscSpaceCreate(PetscObjectComm((PetscObject)sp), &sub);
481: PetscObjectGetName((PetscObject)sp, &name);
482: PetscObjectSetName((PetscObject)sub, name);
483: PetscSpaceSetType(sub, PETSCSPACESUM);
484: PetscSpaceSumSetNumSubspaces(sub, sum->numSumSpaces);
485: PetscSpaceSumSetConcatenate(sub, sum->concatenate);
486: PetscSpaceSetNumComponents(sub, Nc);
487: PetscSpaceSetNumVariables(sub, dim - height);
488: for (PetscInt i = 0; i < sum->numSumSpaces; i++) {
489: PetscSpace subh;
491: PetscSpaceGetHeightSubspace(sum->sumspaces[i], height, &subh);
492: PetscSpaceSumSetSubspace(sub, i, subh);
493: }
494: PetscSpaceSetUp(sub);
495: sum->heightsubspaces[height - 1] = sub;
496: }
497: *subsp = sum->heightsubspaces[height - 1];
498: } else {
499: *subsp = NULL;
500: }
501: return 0;
502: }
504: static PetscErrorCode PetscSpaceInitialize_Sum(PetscSpace sp)
505: {
506: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Sum;
507: sp->ops->setup = PetscSpaceSetUp_Sum;
508: sp->ops->view = PetscSpaceView_Sum;
509: sp->ops->destroy = PetscSpaceDestroy_Sum;
510: sp->ops->getdimension = PetscSpaceGetDimension_Sum;
511: sp->ops->evaluate = PetscSpaceEvaluate_Sum;
512: sp->ops->getheightsubspace = PetscSpaceGetHeightSubspace_Sum;
514: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetNumSubspaces_C", PetscSpaceSumGetNumSubspaces_Sum);
515: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetNumSubspaces_C", PetscSpaceSumSetNumSubspaces_Sum);
516: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetSubspace_C", PetscSpaceSumGetSubspace_Sum);
517: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetSubspace_C", PetscSpaceSumSetSubspace_Sum);
518: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumGetConcatenate_C", PetscSpaceSumGetConcatenate_Sum);
519: PetscObjectComposeFunction((PetscObject)sp, "PetscSpaceSumSetConcatenate_C", PetscSpaceSumSetConcatenate_Sum);
520: return 0;
521: }
523: /*MC
524: PETSCSPACESUM = "sum" - A PetscSpace object that encapsulates a sum of subspaces.
525: That sum can either be direct or concatenate a concatenation.For example if A and B are spaces each with 2 components,
526: 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
527: same number of variables.
529: Level: intermediate
531: .seealso: `PetscSpaceType`, `PetscSpaceCreate()`, `PetscSpaceSetType()`
532: M*/
533: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Sum(PetscSpace sp)
534: {
535: PetscSpace_Sum *sum;
538: PetscNew(&sum);
539: sum->numSumSpaces = PETSC_DEFAULT;
540: sp->data = sum;
541: PetscSpaceInitialize_Sum(sp);
542: return 0;
543: }
545: PETSC_EXTERN PetscErrorCode PetscSpaceCreateSum(PetscInt numSubspaces, const PetscSpace subspaces[], PetscBool concatenate, PetscSpace *sumSpace)
546: {
547: PetscInt i, Nv, Nc = 0;
549: if (sumSpace) PetscSpaceDestroy(sumSpace);
550: PetscSpaceCreate(PetscObjectComm((PetscObject)subspaces[0]), sumSpace);
551: PetscSpaceSetType(*sumSpace, PETSCSPACESUM);
552: PetscSpaceSumSetNumSubspaces(*sumSpace, numSubspaces);
553: PetscSpaceSumSetConcatenate(*sumSpace, concatenate);
554: for (i = 0; i < numSubspaces; ++i) {
555: PetscInt sNc;
557: PetscSpaceSumSetSubspace(*sumSpace, i, subspaces[i]);
558: PetscSpaceGetNumComponents(subspaces[i], &sNc);
559: if (concatenate) Nc += sNc;
560: else Nc = sNc;
561: }
562: PetscSpaceGetNumVariables(subspaces[0], &Nv);
563: PetscSpaceSetNumComponents(*sumSpace, Nc);
564: PetscSpaceSetNumVariables(*sumSpace, Nv);
565: PetscSpaceSetUp(*sumSpace);
567: return 0;
568: }