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, &deg, 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: }