Actual source code: dualspacesum.c

  1: #include <petsc/private/petscfeimpl.h>

  3: /*@
  4:   PetscDualSpaceSumGetNumSubspaces - Get the number of spaces in the sum space

  6:   Input Parameter:
  7: . sp - the dual 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: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumSetNumSubspaces()`
 18: @*/
 19: PetscErrorCode PetscDualSpaceSumGetNumSubspaces(PetscDualSpace sp, PetscInt *numSumSpaces)
 20: {
 21:   PetscFunctionBegin;
 23:   PetscAssertPointer(numSumSpaces, 2);
 24:   PetscTryMethod(sp, "PetscDualSpaceSumGetNumSubspaces_C", (PetscDualSpace, PetscInt *), (sp, numSumSpaces));
 25:   PetscFunctionReturn(PETSC_SUCCESS);
 26: }

 28: /*@
 29:   PetscDualSpaceSumSetNumSubspaces - Set the number of spaces in the sum space

 31:   Input Parameters:
 32: + sp           - the dual 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: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumGetNumSubspaces()`
 41: @*/
 42: PetscErrorCode PetscDualSpaceSumSetNumSubspaces(PetscDualSpace sp, PetscInt numSumSpaces)
 43: {
 44:   PetscFunctionBegin;
 46:   PetscTryMethod(sp, "PetscDualSpaceSumSetNumSubspaces_C", (PetscDualSpace, PetscInt), (sp, numSumSpaces));
 47:   PetscFunctionReturn(PETSC_SUCCESS);
 48: }

 50: /*@
 51:   PetscDualSpaceSumGetConcatenate - Get the concatenate flag for this space.

 53:   Input Parameter:
 54: . sp - the dual space object

 56:   Output Parameter:
 57: . concatenate - flag indicating whether subspaces are concatenated.

 59:   Level: intermediate

 61:   Note:
 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: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumSetConcatenate()`
 67: @*/
 68: PetscErrorCode PetscDualSpaceSumGetConcatenate(PetscDualSpace sp, PetscBool *concatenate)
 69: {
 70:   PetscFunctionBegin;
 72:   PetscTryMethod(sp, "PetscDualSpaceSumGetConcatenate_C", (PetscDualSpace, PetscBool *), (sp, concatenate));
 73:   PetscFunctionReturn(PETSC_SUCCESS);
 74: }

 76: /*@
 77:   PetscDualSpaceSumSetConcatenate - Sets the concatenate flag for this space.

 79:   Input Parameters:
 80: + sp          - the dual 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: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumGetConcatenate()`
 91: @*/
 92: PetscErrorCode PetscDualSpaceSumSetConcatenate(PetscDualSpace sp, PetscBool concatenate)
 93: {
 94:   PetscFunctionBegin;
 96:   PetscTryMethod(sp, "PetscDualSpaceSumSetConcatenate_C", (PetscDualSpace, PetscBool), (sp, concatenate));
 97:   PetscFunctionReturn(PETSC_SUCCESS);
 98: }

100: /*@
101:   PetscDualSpaceSumGetSubspace - Get a space in the sum space

103:   Input Parameters:
104: + sp - the dual space object
105: - s  - The space number

107:   Output Parameter:
108: . subsp - the `PetscDualSpace`

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: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumSetSubspace()`
116: @*/
117: PetscErrorCode PetscDualSpaceSumGetSubspace(PetscDualSpace sp, PetscInt s, PetscDualSpace *subsp)
118: {
119:   PetscFunctionBegin;
121:   PetscAssertPointer(subsp, 3);
122:   PetscTryMethod(sp, "PetscDualSpaceSumGetSubspace_C", (PetscDualSpace, PetscInt, PetscDualSpace *), (sp, s, subsp));
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

126: /*@
127:   PetscDualSpaceSumSetSubspace - Set a space in the sum space

129:   Input Parameters:
130: + sp    - the dual 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: `PETSCDUALSPACESUM`, `PetscDualSpace`, `PetscDualSpaceSumGetSubspace()`
140: @*/
141: PetscErrorCode PetscDualSpaceSumSetSubspace(PetscDualSpace sp, PetscInt s, PetscDualSpace subsp)
142: {
143:   PetscFunctionBegin;
146:   PetscTryMethod(sp, "PetscDualSpaceSumSetSubspace_C", (PetscDualSpace, PetscInt, PetscDualSpace), (sp, s, subsp));
147:   PetscFunctionReturn(PETSC_SUCCESS);
148: }

150: static PetscErrorCode PetscDualSpaceSumGetNumSubspaces_Sum(PetscDualSpace space, PetscInt *numSumSpaces)
151: {
152:   PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)space->data;

154:   PetscFunctionBegin;
155:   *numSumSpaces = sum->numSumSpaces;
156:   PetscFunctionReturn(PETSC_SUCCESS);
157: }

159: static PetscErrorCode PetscDualSpaceSumSetNumSubspaces_Sum(PetscDualSpace space, PetscInt numSumSpaces)
160: {
161:   PetscDualSpace_Sum *sum = (PetscDualSpace_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:   for (PetscInt s = 0; s < Ns; ++s) PetscCall(PetscDualSpaceDestroy(&sum->sumspaces[s]));
168:   PetscCall(PetscFree(sum->sumspaces));

170:   Ns = sum->numSumSpaces = numSumSpaces;
171:   PetscCall(PetscCalloc1(Ns, &sum->sumspaces));
172:   PetscFunctionReturn(PETSC_SUCCESS);
173: }

175: static PetscErrorCode PetscDualSpaceSumGetConcatenate_Sum(PetscDualSpace sp, PetscBool *concatenate)
176: {
177:   PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;

179:   PetscFunctionBegin;
180:   *concatenate = sum->concatenate;
181:   PetscFunctionReturn(PETSC_SUCCESS);
182: }

184: static PetscErrorCode PetscDualSpaceSumSetConcatenate_Sum(PetscDualSpace sp, PetscBool concatenate)
185: {
186:   PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;

188:   PetscFunctionBegin;
189:   PetscCheck(!sum->setupCalled, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Cannot change space concatenation after setup called.");

191:   sum->concatenate = concatenate;
192:   PetscFunctionReturn(PETSC_SUCCESS);
193: }

195: static PetscErrorCode PetscDualSpaceSumGetSubspace_Sum(PetscDualSpace space, PetscInt s, PetscDualSpace *subspace)
196: {
197:   PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)space->data;
198:   PetscInt            Ns  = sum->numSumSpaces;

200:   PetscFunctionBegin;
201:   PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscDualSpaceSumSetNumSubspaces() first");
202:   PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);

204:   *subspace = sum->sumspaces[s];
205:   PetscFunctionReturn(PETSC_SUCCESS);
206: }

208: static PetscErrorCode PetscDualSpaceSumSetSubspace_Sum(PetscDualSpace space, PetscInt s, PetscDualSpace subspace)
209: {
210:   PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)space->data;
211:   PetscInt            Ns  = sum->numSumSpaces;

213:   PetscFunctionBegin;
214:   PetscCheck(!sum->setupCalled, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Cannot change subspace after setup called");
215:   PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_WRONGSTATE, "Must call PetscDualSpaceSumSetNumSubspaces() first");
216:   PetscCheck(s >= 0 && s < Ns, PetscObjectComm((PetscObject)space), PETSC_ERR_ARG_OUTOFRANGE, "Invalid subspace number %" PetscInt_FMT, s);

218:   PetscCall(PetscObjectReference((PetscObject)subspace));
219:   PetscCall(PetscDualSpaceDestroy(&sum->sumspaces[s]));
220:   sum->sumspaces[s] = subspace;
221:   PetscFunctionReturn(PETSC_SUCCESS);
222: }

224: static PetscErrorCode PetscDualSpaceDuplicate_Sum(PetscDualSpace sp, PetscDualSpace spNew)
225: {
226:   PetscInt       num_subspaces, Nc;
227:   PetscBool      concatenate, interleave_basis, interleave_components;
228:   PetscDualSpace subsp_first     = NULL;
229:   PetscDualSpace subsp_dup_first = NULL;
230:   DM             K;

232:   PetscFunctionBegin;
233:   PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &num_subspaces));
234:   PetscCall(PetscDualSpaceSumSetNumSubspaces(spNew, num_subspaces));
235:   PetscCall(PetscDualSpaceSumGetConcatenate(sp, &concatenate));
236:   PetscCall(PetscDualSpaceSumSetConcatenate(spNew, concatenate));
237:   PetscCall(PetscDualSpaceSumGetInterleave(sp, &interleave_basis, &interleave_components));
238:   PetscCall(PetscDualSpaceSumSetInterleave(spNew, interleave_basis, interleave_components));
239:   PetscCall(PetscDualSpaceGetDM(sp, &K));
240:   PetscCall(PetscDualSpaceSetDM(spNew, K));
241:   PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
242:   PetscCall(PetscDualSpaceSetNumComponents(spNew, Nc));
243:   for (PetscInt s = 0; s < num_subspaces; s++) {
244:     PetscDualSpace subsp, subspNew;

246:     PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
247:     if (s == 0) {
248:       subsp_first = subsp;
249:       PetscCall(PetscDualSpaceDuplicate(subsp, &subsp_dup_first));
250:       subspNew = subsp_dup_first;
251:     } else if (subsp == subsp_first) {
252:       PetscCall(PetscObjectReference((PetscObject)subsp_dup_first));
253:       subspNew = subsp_dup_first;
254:     } else {
255:       PetscCall(PetscDualSpaceDuplicate(subsp, &subspNew));
256:     }
257:     PetscCall(PetscDualSpaceSumSetSubspace(spNew, s, subspNew));
258:     PetscCall(PetscDualSpaceDestroy(&subspNew));
259:   }
260:   PetscFunctionReturn(PETSC_SUCCESS);
261: }

263: static PetscErrorCode PetscDualSpaceSumCreateQuadrature(PetscInt Ns, PetscInt cdim, PetscBool uniform_points, PetscQuadrature subquads[], PetscQuadrature *fullquad)
264: {
265:   PetscReal *points;
266:   PetscInt   Npoints;

268:   PetscFunctionBegin;
269:   if (uniform_points) {
270:     PetscCall(PetscObjectReference((PetscObject)subquads[0]));
271:     *fullquad = subquads[0];
272:     PetscFunctionReturn(PETSC_SUCCESS);
273:   }
274:   Npoints = 0;
275:   for (PetscInt s = 0; s < Ns; s++) {
276:     PetscInt subNpoints;

278:     if (!subquads[s]) continue;
279:     PetscCall(PetscQuadratureGetData(subquads[s], NULL, NULL, &subNpoints, NULL, NULL));
280:     Npoints += subNpoints;
281:   }
282:   *fullquad = NULL;
283:   if (!Npoints) PetscFunctionReturn(PETSC_SUCCESS);
284:   PetscCall(PetscMalloc1(Npoints * cdim, &points));
285:   for (PetscInt s = 0, offset = 0; s < Ns; s++) {
286:     PetscInt         subNpoints;
287:     const PetscReal *subpoints;

289:     PetscCall(PetscQuadratureGetData(subquads[s], NULL, NULL, &subNpoints, &subpoints, NULL));
290:     PetscCall(PetscArraycpy(&points[offset], subpoints, cdim * subNpoints));
291:     offset += cdim * subNpoints;
292:   }
293:   PetscCall(PetscQuadratureCreate(PETSC_COMM_SELF, fullquad));
294:   PetscCall(PetscQuadratureSetData(*fullquad, cdim, 0, Npoints, points, NULL));
295:   PetscFunctionReturn(PETSC_SUCCESS);
296: }

298: static PetscErrorCode PetscDualSpaceSumCreateMatrix(PetscDualSpace sp, Mat submats[], ISLocalToGlobalMapping map_rows[], ISLocalToGlobalMapping map_cols[], PetscQuadrature fullquad, Mat *fullmat)
299: {
300:   Mat          mat;
301:   PetscInt    *i = NULL, *j = NULL;
302:   PetscScalar *v = NULL;
303:   PetscInt     npoints;
304:   PetscInt     nrows = 0, ncols, nnz = 0;
305:   PetscInt     Nc, num_subspaces;

307:   PetscFunctionBegin;
308:   PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
309:   PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &num_subspaces));
310:   PetscCall(PetscQuadratureGetData(fullquad, NULL, NULL, &npoints, NULL, NULL));
311:   ncols = Nc * npoints;
312:   for (PetscInt s = 0; s < num_subspaces; s++) {
313:     // Get the COO data for each matrix, map the is and js, and append to growing COO data
314:     PetscInt        sNrows, sNcols;
315:     Mat             smat;
316:     const PetscInt *si;
317:     const PetscInt *sj;
318:     PetscScalar    *sv;
319:     PetscMemType    memtype;
320:     PetscInt        snz;
321:     PetscInt        snz_actual;
322:     PetscInt       *cooi;
323:     PetscInt       *cooj;
324:     PetscScalar    *coov;
325:     PetscScalar    *v_new;
326:     PetscInt       *i_new;
327:     PetscInt       *j_new;

329:     if (!submats[s]) continue;
330:     PetscCall(MatGetSize(submats[s], &sNrows, &sNcols));
331:     nrows += sNrows;
332:     PetscCall(MatConvert(submats[s], MATSEQAIJ, MAT_INITIAL_MATRIX, &smat));
333:     PetscCall(MatBindToCPU(smat, PETSC_TRUE));
334:     PetscCall(MatSeqAIJGetCSRAndMemType(smat, &si, &sj, &sv, &memtype));
335:     PetscCheck(memtype == PETSC_MEMTYPE_HOST, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Could not convert matrix to host memory");
336:     snz = si[sNrows];

338:     PetscCall(PetscMalloc1(nnz + snz, &v_new));
339:     PetscCall(PetscArraycpy(v_new, v, nnz));
340:     PetscCall(PetscFree(v));
341:     v = v_new;

343:     PetscCall(PetscMalloc1(nnz + snz, &i_new));
344:     PetscCall(PetscArraycpy(i_new, i, nnz));
345:     PetscCall(PetscFree(i));
346:     i = i_new;

348:     PetscCall(PetscMalloc1(nnz + snz, &j_new));
349:     PetscCall(PetscArraycpy(j_new, j, nnz));
350:     PetscCall(PetscFree(j));
351:     j = j_new;

353:     PetscCall(PetscMalloc2(snz, &cooi, snz, &cooj));

355:     coov = &v[nnz];

357:     snz_actual = 0;
358:     for (PetscInt row = 0; row < sNrows; row++) {
359:       for (PetscInt k = si[row]; k < si[row + 1]; k++, snz_actual++) {
360:         cooi[snz_actual] = row;
361:         cooj[snz_actual] = sj[k];
362:         coov[snz_actual] = sv[k];
363:       }
364:     }
365:     PetscCall(MatDestroy(&smat));

367:     if (snz_actual > 0) {
368:       PetscCall(ISLocalToGlobalMappingApply(map_rows[s], snz_actual, cooi, &i[nnz]));
369:       PetscCall(ISLocalToGlobalMappingApply(map_cols[s], snz_actual, cooj, &j[nnz]));
370:       nnz += snz_actual;
371:     }
372:     PetscCall(PetscFree2(cooi, cooj));
373:   }
374:   PetscCall(MatCreate(PETSC_COMM_SELF, fullmat));
375:   mat = *fullmat;
376:   PetscCall(MatSetSizes(mat, nrows, ncols, nrows, ncols));
377:   PetscCall(MatSetType(mat, MATSEQAIJ));
378:   PetscCall(MatSetPreallocationCOO(mat, nnz, i, j));
379:   PetscCall(MatSetValuesCOO(mat, v, INSERT_VALUES));
380:   PetscCall(PetscFree(i));
381:   PetscCall(PetscFree(j));
382:   PetscCall(PetscFree(v));
383:   PetscFunctionReturn(PETSC_SUCCESS);
384: }

386: static PetscErrorCode PetscDualSpaceSumCreateMappings(PetscDualSpace sp, PetscBool interior, PetscBool uniform_points, ISLocalToGlobalMapping map_row[], ISLocalToGlobalMapping map_col[])
387: {
388:   PetscSection section;
389:   PetscInt     pStart, pEnd, Ns, Nc;
390:   PetscInt     num_points = 0;
391:   PetscBool    interleave_basis, interleave_components, concatenate;
392:   PetscInt    *roffset;

394:   PetscFunctionBegin;
395:   if (interior) {
396:     PetscCall(PetscDualSpaceGetInteriorSection(sp, &section));
397:   } else {
398:     PetscCall(PetscDualSpaceGetSection(sp, &section));
399:   }
400:   PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
401:   PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
402:   PetscCall(PetscDualSpaceSumGetInterleave(sp, &interleave_basis, &interleave_components));
403:   PetscCall(PetscDualSpaceSumGetConcatenate(sp, &concatenate));
404:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
405:   for (PetscInt p = pStart; p < pEnd; p++) {
406:     PetscInt dof;

408:     PetscCall(PetscSectionGetDof(section, p, &dof));
409:     num_points += (dof > 0);
410:   }
411:   PetscCall(PetscCalloc1(pEnd - pStart, &roffset));
412:   for (PetscInt s = 0, coffset = 0; s < Ns; s++) {
413:     PetscDualSpace  subsp;
414:     PetscSection    subsection;
415:     IS              is_row, is_col;
416:     PetscInt        subNb;
417:     PetscInt        sNc, sNpoints, sNcols;
418:     PetscQuadrature quad;

420:     PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
421:     PetscCall(PetscDualSpaceGetNumComponents(subsp, &sNc));
422:     if (interior) {
423:       PetscCall(PetscDualSpaceGetInteriorSection(subsp, &subsection));
424:       PetscCall(PetscDualSpaceGetInteriorData(subsp, &quad, NULL));
425:     } else {
426:       PetscCall(PetscDualSpaceGetSection(subsp, &subsection));
427:       PetscCall(PetscDualSpaceGetAllData(subsp, &quad, NULL));
428:     }
429:     PetscCall(PetscSectionGetStorageSize(subsection, &subNb));
430:     if (num_points == 1) {
431:       PetscInt rstride;

433:       rstride = interleave_basis ? Ns : 1;
434:       PetscCall(ISCreateStride(PETSC_COMM_SELF, subNb, roffset[0], rstride, &is_row));
435:       roffset[0] += interleave_basis ? 1 : subNb;
436:     } else {
437:       PetscInt *rows;

439:       PetscCall(PetscMalloc1(subNb, &rows));
440:       for (PetscInt p = pStart; p < pEnd; p++) {
441:         PetscInt subdof, dof, off, suboff, stride;

443:         PetscCall(PetscSectionGetOffset(subsection, p, &suboff));
444:         PetscCall(PetscSectionGetDof(subsection, p, &subdof));
445:         PetscCall(PetscSectionGetOffset(section, p, &off));
446:         PetscCall(PetscSectionGetDof(section, p, &dof));
447:         PetscCheck(subdof * Ns == dof || !interleave_basis, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Basis cannot be interleaved");
448:         stride = interleave_basis ? Ns : 1;
449:         for (PetscInt k = 0; k < subdof; k++) { rows[suboff + k] = off + roffset[p - pStart] + k * stride; }
450:         roffset[p - pStart] += interleave_basis ? 1 : subdof;
451:       }
452:       PetscCall(ISCreateGeneral(PETSC_COMM_SELF, subNb, rows, PETSC_OWN_POINTER, &is_row));
453:     }

455:     sNpoints = 0;
456:     if (quad) PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &sNpoints, NULL, NULL));
457:     sNcols = sNpoints * sNc;

459:     if (!concatenate) {
460:       PetscCall(ISCreateStride(PETSC_COMM_SELF, sNcols, coffset, 1, &is_col));
461:       coffset += uniform_points ? 0 : sNcols;
462:     } else {
463:       if (uniform_points && interleave_components) {
464:         PetscCall(ISCreateStride(PETSC_COMM_SELF, sNcols, coffset, Ns, &is_col));
465:         coffset += 1;
466:       } else {
467:         PetscInt *cols;

469:         PetscCall(PetscMalloc1(sNcols, &cols));
470:         for (PetscInt p = 0, r = 0; p < sNpoints; p++) {
471:           for (PetscInt c = 0; c < sNc; c++, r++) { cols[r] = coffset + p * Nc + c; }
472:         }
473:         coffset += uniform_points ? sNc : Nc * sNpoints + sNc;
474:         PetscCall(ISCreateGeneral(PETSC_COMM_SELF, sNcols, cols, PETSC_OWN_POINTER, &is_col));
475:       }
476:     }
477:     PetscCall(ISLocalToGlobalMappingCreateIS(is_row, &map_row[s]));
478:     PetscCall(ISLocalToGlobalMappingCreateIS(is_col, &map_col[s]));
479:     PetscCall(ISDestroy(&is_row));
480:     PetscCall(ISDestroy(&is_col));
481:   }
482:   PetscCall(PetscFree(roffset));
483:   PetscFunctionReturn(PETSC_SUCCESS);
484: }

486: static PetscErrorCode PetscDualSpaceCreateFacetSubspace_Sum(PetscDualSpace sp, PetscInt f, PetscDualSpace *bdsp)
487: {
488:   PetscInt       k, Nc, Nk, Nknew, Ncnew, Ns;
489:   PetscInt       dim, pointDim = -1;
490:   PetscInt       depth, Ncopies;
491:   PetscBool      any;
492:   DM             dm, K;
493:   DMPolytopeType ct;

495:   PetscFunctionBegin;
496:   PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
497:   any = PETSC_FALSE;
498:   for (PetscInt s = 0; s < Ns; s++) {
499:     PetscDualSpace subsp, subpointsp;

501:     PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
502:     PetscCall(PetscDualSpaceGetPointSubspace(subsp, f, &subpointsp));
503:     if (subpointsp) any = PETSC_TRUE;
504:   }
505:   if (!any) {
506:     *bdsp = NULL;
507:     PetscFunctionReturn(PETSC_SUCCESS);
508:   }
509:   PetscCall(PetscDualSpaceGetDM(sp, &dm));
510:   PetscCall(DMGetDimension(dm, &dim));
511:   PetscCall(DMPlexGetDepth(dm, &depth));
512:   PetscCall(PetscDualSpaceDuplicate(sp, bdsp));
513:   PetscCheck((depth == dim) || (depth == 1), PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unsupported interpolation state of reference element");
514:   pointDim = (depth == dim) ? (dim - 1) : 0;
515:   PetscCall(DMPlexGetCellType(dm, f, &ct));
516:   PetscCall(DMPlexCreateReferenceCell(PETSC_COMM_SELF, ct, &K));
517:   PetscCall(PetscDualSpaceSetDM(*bdsp, K));
518:   PetscCall(PetscDualSpaceGetFormDegree(sp, &k));
519:   PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
520:   PetscCall(PetscDTBinomialInt(dim, PetscAbsInt(k), &Nk));
521:   Ncopies = Nc / Nk;
522:   PetscCall(PetscDTBinomialInt(pointDim, PetscAbsInt(k), &Nknew));
523:   Ncnew = Nknew * Ncopies;
524:   PetscCall(PetscDualSpaceSetNumComponents(*bdsp, Ncnew));
525:   for (PetscInt s = 0; s < Ns; s++) {
526:     PetscDualSpace subsp, subpointsp;

528:     PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
529:     PetscCall(PetscDualSpaceGetPointSubspace(subsp, f, &subpointsp));
530:     if (subpointsp) {
531:       PetscCall(PetscObjectReference((PetscObject)subpointsp));
532:     } else {
533:       // make an empty dual space
534:       PetscInt subNc, subNcopies, subpointNc;

536:       PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)subsp), &subpointsp));
537:       PetscCall(PetscDualSpaceGetNumComponents(subsp, &subNc));
538:       subNcopies = subNc / Nk;
539:       subpointNc = subNcopies * Nknew;
540:       PetscCall(PetscDualSpaceSetType(subpointsp, PETSCDUALSPACESIMPLE));
541:       PetscCall(PetscDualSpaceSimpleSetDimension(subpointsp, 0));
542:       PetscCall(PetscDualSpaceSetFormDegree(subpointsp, k));
543:       PetscCall(PetscDualSpaceSetNumComponents(subpointsp, subpointNc));
544:     }
545:     // K should be equal to subpointsp->dm, but we want it to be exactly the same
546:     PetscCall(PetscObjectReference((PetscObject)K));
547:     PetscCall(DMDestroy(&subpointsp->dm));
548:     subpointsp->dm = K;
549:     PetscCall(PetscDualSpaceSetUp(subpointsp));
550:     PetscCall(PetscDualSpaceSumSetSubspace(*bdsp, s, subpointsp));
551:     PetscCall(PetscDualSpaceDestroy(&subpointsp));
552:   }
553:   PetscCall(DMDestroy(&K));
554:   PetscCall(PetscDualSpaceSetUp(*bdsp));
555:   PetscFunctionReturn(PETSC_SUCCESS);
556: }

558: static PetscErrorCode PetscDualSpaceSumIsUniform(PetscDualSpace sp, PetscBool *is_uniform)
559: {
560:   PetscDualSpace_Sum *sum     = (PetscDualSpace_Sum *)sp->data;
561:   PetscBool           uniform = PETSC_TRUE;

563:   PetscFunctionBegin;
564:   for (PetscInt s = 1; s < sum->numSumSpaces; s++) {
565:     if (sum->sumspaces[s] != sum->sumspaces[0]) {
566:       uniform = PETSC_FALSE;
567:       break;
568:     }
569:   }
570:   *is_uniform = uniform;
571:   PetscFunctionReturn(PETSC_SUCCESS);
572: }

574: static PetscErrorCode PetscDualSpaceGetSymmetries_Sum(PetscDualSpace sp, const PetscInt ****perms, const PetscScalar ****flips)
575: {
576:   PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;

578:   PetscFunctionBegin;
579:   if (!sum->symComputed) {
580:     PetscInt       Ns;
581:     PetscBool      any_perms = PETSC_FALSE;
582:     PetscBool      any_flips = PETSC_FALSE;
583:     PetscInt    ***symperms  = NULL;
584:     PetscScalar ***symflips  = NULL;

586:     sum->symComputed = PETSC_TRUE;
587:     PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
588:     for (PetscInt s = 0; s < Ns; s++) {
589:       PetscDualSpace       subsp;
590:       const PetscInt    ***sub_perms;
591:       const PetscScalar ***sub_flips;

593:       PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
594:       PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
595:       if (sub_perms) any_perms = PETSC_TRUE;
596:       if (sub_flips) any_flips = PETSC_TRUE;
597:     }
598:     if (any_perms || any_flips) {
599:       DM       K;
600:       PetscInt pStart, pEnd, numPoints;
601:       PetscInt spintdim;

603:       PetscCall(PetscDualSpaceGetDM(sp, &K));
604:       PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
605:       numPoints = pEnd - pStart;
606:       PetscCall(PetscCalloc1(numPoints, &symperms));
607:       PetscCall(PetscCalloc1(numPoints, &symflips));
608:       PetscCall(PetscDualSpaceGetBoundarySymmetries_Internal(sp, symperms, symflips));
609:       // get interior symmetries
610:       PetscCall(PetscDualSpaceGetInteriorDimension(sp, &spintdim));
611:       if (spintdim) {
612:         PetscInt       groupSize;
613:         PetscInt     **cellPerms;
614:         PetscScalar  **cellFlips;
615:         DMPolytopeType ct;

617:         PetscCall(DMPlexGetCellType(K, 0, &ct));
618:         groupSize       = DMPolytopeTypeGetNumArrangements(ct);
619:         sum->numSelfSym = groupSize;
620:         sum->selfSymOff = groupSize / 2;
621:         PetscCall(PetscCalloc1(groupSize, &cellPerms));
622:         PetscCall(PetscCalloc1(groupSize, &cellFlips));
623:         symperms[0] = &cellPerms[groupSize / 2];
624:         symflips[0] = &cellFlips[groupSize / 2];
625:         for (PetscInt o = -groupSize / 2; o < groupSize / 2; o++) {
626:           PetscBool any_o_perms = PETSC_FALSE;
627:           PetscBool any_o_flips = PETSC_FALSE;

629:           for (PetscInt s = 0; s < Ns; s++) {
630:             PetscDualSpace       subsp;
631:             const PetscInt    ***sub_perms;
632:             const PetscScalar ***sub_flips;

634:             PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
635:             PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
636:             if (sub_perms && sub_perms[0] && sub_perms[0][o]) any_o_perms = PETSC_TRUE;
637:             if (sub_flips && sub_flips[0] && sub_flips[0][o]) any_o_flips = PETSC_TRUE;
638:           }
639:           if (any_o_perms) {
640:             PetscInt *o_perm;

642:             PetscCall(PetscMalloc1(spintdim, &o_perm));
643:             for (PetscInt i = 0; i < spintdim; i++) o_perm[i] = i;
644:             for (PetscInt s = 0; s < Ns; s++) {
645:               PetscDualSpace       subsp;
646:               const PetscInt    ***sub_perms;
647:               const PetscScalar ***sub_flips;

649:               PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
650:               PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
651:               if (sub_perms && sub_perms[0] && sub_perms[0][o]) {
652:                 PetscInt  subspdim;
653:                 PetscInt *range, *domain;
654:                 PetscInt *range_mapped, *domain_mapped;

656:                 PetscCall(PetscDualSpaceGetInteriorDimension(subsp, &subspdim));
657:                 PetscCall(PetscMalloc4(subspdim, &range, subspdim, &range_mapped, subspdim, &domain, subspdim, &domain_mapped));
658:                 for (PetscInt i = 0; i < subspdim; i++) domain[i] = i;
659:                 PetscCall(PetscArraycpy(range, sub_perms[0][o], subspdim));
660:                 PetscCall(ISLocalToGlobalMappingApply(sum->int_rows[s], subspdim, domain, domain_mapped));
661:                 PetscCall(ISLocalToGlobalMappingApply(sum->int_rows[s], subspdim, range, range_mapped));
662:                 for (PetscInt i = 0; i < subspdim; i++) o_perm[domain_mapped[i]] = range_mapped[i];
663:                 PetscCall(PetscFree4(range, range_mapped, domain, domain_mapped));
664:               }
665:             }
666:             symperms[0][o] = o_perm;
667:           }
668:           if (any_o_flips) {
669:             PetscScalar *o_flip;

671:             PetscCall(PetscMalloc1(spintdim, &o_flip));
672:             for (PetscInt i = 0; i < spintdim; i++) o_flip[i] = 1.0;
673:             for (PetscInt s = 0; s < Ns; s++) {
674:               PetscDualSpace       subsp;
675:               const PetscInt    ***sub_perms;
676:               const PetscScalar ***sub_flips;

678:               PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
679:               PetscCall(PetscDualSpaceGetSymmetries(subsp, &sub_perms, &sub_flips));
680:               if (sub_perms && sub_perms[0] && sub_perms[0][o]) {
681:                 PetscInt  subspdim;
682:                 PetscInt *domain;
683:                 PetscInt *domain_mapped;

685:                 PetscCall(PetscDualSpaceGetInteriorDimension(subsp, &subspdim));
686:                 PetscCall(PetscMalloc2(subspdim, &domain, subspdim, &domain_mapped));
687:                 for (PetscInt i = 0; i < subspdim; i++) domain[i] = i;
688:                 PetscCall(ISLocalToGlobalMappingApply(sum->int_rows[s], subspdim, domain, domain_mapped));
689:                 for (PetscInt i = 0; i < subspdim; i++) o_flip[domain_mapped[i]] = sub_perms[0][o][i];
690:                 PetscCall(PetscFree2(domain, domain_mapped));
691:               }
692:             }
693:             symflips[0][o] = o_flip;
694:           }
695:         }
696:         {
697:           PetscBool any_perms = PETSC_FALSE;
698:           PetscBool any_flips = PETSC_FALSE;
699:           for (PetscInt o = -groupSize / 2; o < groupSize / 2; o++) {
700:             if (symperms[0][o]) any_perms = PETSC_TRUE;
701:             if (symflips[0][o]) any_flips = PETSC_TRUE;
702:           }
703:           if (!any_perms) {
704:             PetscCall(PetscFree(cellPerms));
705:             symperms[0] = NULL;
706:           }
707:           if (!any_flips) {
708:             PetscCall(PetscFree(cellFlips));
709:             symflips[0] = NULL;
710:           }
711:         }
712:       }
713:       if (!any_perms) {
714:         PetscCall(PetscFree(symperms));
715:         symperms = NULL;
716:       }
717:       if (!any_flips) {
718:         PetscCall(PetscFree(symflips));
719:         symflips = NULL;
720:       }
721:     }
722:     sum->symperms = symperms;
723:     sum->symflips = symflips;
724:   }
725:   if (perms) *perms = (const PetscInt ***)sum->symperms;
726:   if (flips) *flips = (const PetscScalar ***)sum->symflips;
727:   PetscFunctionReturn(PETSC_SUCCESS);
728: }

730: static PetscErrorCode PetscDualSpaceSetUp_Sum(PetscDualSpace sp)
731: {
732:   PetscDualSpace_Sum *sum         = (PetscDualSpace_Sum *)sp->data;
733:   PetscBool           concatenate = PETSC_TRUE;
734:   PetscBool           uniform;
735:   PetscInt            Ns, Nc, i, sum_Nc = 0;
736:   PetscInt            minNc, maxNc;
737:   PetscInt            minForm, maxForm, cdim, depth;
738:   DM                  K;
739:   PetscQuadrature    *all_quads = NULL;
740:   PetscQuadrature    *int_quads = NULL;
741:   Mat                *all_mats  = NULL;
742:   Mat                *int_mats  = NULL;

744:   PetscFunctionBegin;
745:   if (sum->setupCalled) PetscFunctionReturn(PETSC_SUCCESS);
746:   sum->setupCalled = PETSC_TRUE;

748:   PetscCall(PetscDualSpaceSumGetNumSubspaces(sp, &Ns));
749:   PetscCheck(Ns >= 0, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Cannot have %" PetscInt_FMT " subspaces", Ns);

751:   // step 1: make sure they share a DM
752:   PetscCall(PetscDualSpaceGetDM(sp, &K));
753:   PetscCall(DMGetCoordinateDim(K, &cdim));
754:   PetscCall(DMPlexGetDepth(K, &depth));
755:   PetscCall(PetscDualSpaceSumIsUniform(sp, &sp->uniform));
756:   uniform = sp->uniform;
757:   {
758:     for (PetscInt s = 0; s < Ns; s++) {
759:       PetscDualSpace subsp;
760:       DM             sub_K;

762:       PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
763:       PetscCall(PetscDualSpaceSetUp(subsp));
764:       PetscCall(PetscDualSpaceGetDM(subsp, &sub_K));
765:       if (s == 0 && K == NULL) {
766:         PetscCall(PetscDualSpaceSetDM(sp, sub_K));
767:         K = sub_K;
768:       }
769:       PetscCheck(sub_K == K, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Subspace %" PetscInt_FMT " does not have the same DM as the sum space", s);
770:     }
771:   }

773:   // step 2: count components
774:   PetscCall(PetscDualSpaceGetNumComponents(sp, &Nc));
775:   PetscCall(PetscDualSpaceSumGetConcatenate(sp, &concatenate));
776:   minNc   = Nc;
777:   maxNc   = Nc;
778:   minForm = sp->k == PETSC_FORM_DEGREE_UNDEFINED ? PETSC_INT_MAX : sp->k;
779:   maxForm = sp->k == PETSC_FORM_DEGREE_UNDEFINED ? PETSC_INT_MIN : sp->k;
780:   for (i = 0; i < Ns; ++i) {
781:     PetscInt       sNc, formDegree;
782:     PetscDualSpace si;

784:     PetscCall(PetscDualSpaceSumGetSubspace(sp, i, &si));
785:     PetscCall(PetscDualSpaceSetUp(si));
786:     PetscCall(PetscDualSpaceGetNumComponents(si, &sNc));
787:     if (sNc != Nc) PetscCheck(concatenate, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_WRONGSTATE, "Subspace as a different number of components but space does not concatenate components");
788:     minNc = PetscMin(minNc, sNc);
789:     maxNc = PetscMax(maxNc, sNc);
790:     sum_Nc += sNc;
791:     PetscCall(PetscDualSpaceGetFormDegree(si, &formDegree));
792:     minForm = PetscMin(minForm, formDegree);
793:     maxForm = PetscMax(maxForm, formDegree);
794:   }
795:   sp->k = (minForm != maxForm) ? PETSC_FORM_DEGREE_UNDEFINED : minForm;

797:   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);
798:   else PetscCheck(minNc == Nc && maxNc == Nc, PetscObjectComm((PetscObject)sp), PETSC_ERR_ARG_OUTOFRANGE, "Subspaces must have same number of components as the target space.");

800:   /* A PetscDualSpace should have a fixed number of components, but
801:      if the spaces we are combining have different form degrees, they will not
802:      have the same number of components on subcomponents of the boundary,
803:      so we do not try to create boundary dual spaces in this case */
804:   if (sp->k != PETSC_FORM_DEGREE_UNDEFINED && depth > 0) {
805:     PetscInt  pStart, pEnd;
806:     PetscInt *pStratStart, *pStratEnd;

808:     PetscCall(DMPlexGetDepth(K, &depth));
809:     PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
810:     PetscCall(PetscCalloc1(pEnd, &sp->pointSpaces));
811:     PetscCall(PetscMalloc2(depth + 1, &pStratStart, depth + 1, &pStratEnd));
812:     for (PetscInt d = 0; d <= depth; ++d) PetscCall(DMPlexGetDepthStratum(K, d, &pStratStart[d], &pStratEnd[d]));

814:     for (PetscInt p = pStratStart[depth - 1]; p < pStratEnd[depth - 1]; p++) { /* calculate the facet dual spaces */
815:       PetscReal      v0[3];
816:       DMPolytopeType ptype;
817:       PetscReal      J[9], detJ;
818:       PetscInt       q;

820:       PetscCall(DMPlexComputeCellGeometryAffineFEM(K, p, v0, J, NULL, &detJ));
821:       PetscCall(DMPlexGetCellType(K, p, &ptype));

823:       /* compare to previous facets: if computed, reference that dualspace */
824:       for (q = pStratStart[depth - 1]; q < p; q++) {
825:         DMPolytopeType qtype;

827:         PetscCall(DMPlexGetCellType(K, q, &qtype));
828:         if (qtype == ptype) break;
829:       }
830:       if (q < p) { /* this facet has the same dual space as that one */
831:         PetscCall(PetscObjectReference((PetscObject)sp->pointSpaces[q]));
832:         sp->pointSpaces[p] = sp->pointSpaces[q];
833:         continue;
834:       }
835:       /* if not, recursively compute this dual space */
836:       PetscCall(PetscDualSpaceCreateFacetSubspace_Sum(sp, p, &sp->pointSpaces[p]));
837:     }
838:     for (PetscInt h = 2; h <= depth; h++) { /* get the higher subspaces from the facet subspaces */
839:       PetscInt hd = depth - h;

841:       for (PetscInt p = pStratStart[hd]; p < pStratEnd[hd]; p++) {
842:         PetscInt        suppSize;
843:         const PetscInt *supp;
844:         DM              qdm;
845:         PetscDualSpace  qsp, psp;
846:         PetscInt        c, coneSize, q;
847:         const PetscInt *cone;
848:         const PetscInt *refCone;

850:         PetscCall(DMPlexGetSupportSize(K, p, &suppSize));
851:         PetscCall(DMPlexGetSupport(K, p, &supp));
852:         q   = supp[0];
853:         qsp = sp->pointSpaces[q];
854:         PetscCall(DMPlexGetConeSize(K, q, &coneSize));
855:         PetscCall(DMPlexGetCone(K, q, &cone));
856:         for (c = 0; c < coneSize; c++)
857:           if (cone[c] == p) break;
858:         PetscCheck(c != coneSize, PetscObjectComm((PetscObject)K), PETSC_ERR_PLIB, "cone/support mismatch");
859:         if (!qsp) {
860:           sp->pointSpaces[p] = NULL;
861:           continue;
862:         }
863:         PetscCall(PetscDualSpaceGetDM(qsp, &qdm));
864:         PetscCall(DMPlexGetCone(qdm, 0, &refCone));
865:         /* get the equivalent dual space from the support dual space */
866:         PetscCall(PetscDualSpaceGetPointSubspace(qsp, refCone[c], &psp));
867:         PetscCall(PetscObjectReference((PetscObject)psp));
868:         sp->pointSpaces[p] = psp;
869:       }
870:     }
871:     PetscCall(PetscFree2(pStratStart, pStratEnd));
872:   }

874:   sum->uniform = uniform;
875:   PetscCall(PetscCalloc1(Ns, &sum->all_rows));
876:   PetscCall(PetscCalloc1(Ns, &sum->all_cols));
877:   PetscCall(PetscCalloc1(Ns, &sum->int_rows));
878:   PetscCall(PetscCalloc1(Ns, &sum->int_cols));
879:   PetscCall(PetscMalloc4(Ns, &all_quads, Ns, &all_mats, Ns, &int_quads, Ns, &int_mats));
880:   {
881:     // test for uniform all points and uniform interior points
882:     PetscBool       uniform_all         = PETSC_TRUE;
883:     PetscBool       uniform_interior    = PETSC_TRUE;
884:     PetscQuadrature quad_all_first      = NULL;
885:     PetscQuadrature quad_interior_first = NULL;
886:     PetscInt        pStart, pEnd;

888:     PetscCall(DMPlexGetChart(K, &pStart, &pEnd));
889:     PetscCall(PetscDualSpaceSectionCreate_Internal(sp, &sp->pointSection));

891:     for (PetscInt p = pStart; p < pEnd; p++) {
892:       PetscInt full_dof = 0;

894:       for (PetscInt s = 0; s < Ns; s++) {
895:         PetscDualSpace subsp;
896:         PetscSection   subsection;
897:         PetscInt       dof;

899:         PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
900:         PetscCall(PetscDualSpaceGetSection(subsp, &subsection));
901:         PetscCall(PetscSectionGetDof(subsection, p, &dof));
902:         full_dof += dof;
903:       }
904:       PetscCall(PetscSectionSetDof(sp->pointSection, p, full_dof));
905:     }
906:     PetscCall(PetscDualSpaceSectionSetUp_Internal(sp, sp->pointSection));

908:     for (PetscInt s = 0; s < Ns; s++) {
909:       PetscDualSpace  subsp;
910:       PetscQuadrature subquad_all;
911:       PetscQuadrature subquad_interior;
912:       Mat             submat_all;
913:       Mat             submat_interior;

915:       PetscCall(PetscDualSpaceSumGetSubspace(sp, s, &subsp));
916:       PetscCall(PetscDualSpaceGetAllData(subsp, &subquad_all, &submat_all));
917:       PetscCall(PetscDualSpaceGetInteriorData(subsp, &subquad_interior, &submat_interior));
918:       if (!s) {
919:         quad_all_first      = subquad_all;
920:         quad_interior_first = subquad_interior;
921:       } else {
922:         if (subquad_all != quad_all_first) uniform_all = PETSC_FALSE;
923:         if (subquad_interior != quad_interior_first) uniform_interior = PETSC_FALSE;
924:       }
925:       all_quads[s] = subquad_all;
926:       int_quads[s] = subquad_interior;
927:       all_mats[s]  = submat_all;
928:       int_mats[s]  = submat_interior;
929:     }
930:     sum->uniform_all_points      = uniform_all;
931:     sum->uniform_interior_points = uniform_interior;

933:     PetscCall(PetscDualSpaceSumCreateMappings(sp, PETSC_TRUE, uniform_interior, sum->int_rows, sum->int_cols));
934:     PetscCall(PetscDualSpaceSumCreateQuadrature(Ns, cdim, uniform_interior, int_quads, &sp->intNodes));
935:     if (sp->intNodes) PetscCall(PetscDualSpaceSumCreateMatrix(sp, int_mats, sum->int_rows, sum->int_cols, sp->intNodes, &sp->intMat));

937:     PetscCall(PetscDualSpaceSumCreateMappings(sp, PETSC_FALSE, uniform_all, sum->all_rows, sum->all_cols));
938:     PetscCall(PetscDualSpaceSumCreateQuadrature(Ns, cdim, uniform_all, all_quads, &sp->allNodes));
939:     if (sp->allNodes) PetscCall(PetscDualSpaceSumCreateMatrix(sp, all_mats, sum->all_rows, sum->all_cols, sp->allNodes, &sp->allMat));
940:   }
941:   PetscCall(PetscFree4(all_quads, all_mats, int_quads, int_mats));
942:   PetscCall(PetscDualSpaceComputeFunctionalsFromAllData(sp));
943:   PetscFunctionReturn(PETSC_SUCCESS);
944: }

946: static PetscErrorCode PetscDualSpaceSumView_Ascii(PetscDualSpace sp, PetscViewer v)
947: {
948:   PetscDualSpace_Sum *sum         = (PetscDualSpace_Sum *)sp->data;
949:   PetscBool           concatenate = sum->concatenate;
950:   PetscInt            i, Ns = sum->numSumSpaces;

952:   PetscFunctionBegin;
953:   if (concatenate) PetscCall(PetscViewerASCIIPrintf(v, "Sum dual space of %" PetscInt_FMT " concatenated subspaces%s\n", Ns, sum->uniform ? " (all identical)" : ""));
954:   else PetscCall(PetscViewerASCIIPrintf(v, "Sum dual space of %" PetscInt_FMT " subspaces%s\n", Ns, sum->uniform ? " (all identical)" : ""));
955:   for (i = 0; i < (sum->uniform ? (Ns > 0 ? 1 : 0) : Ns); ++i) {
956:     PetscCall(PetscViewerASCIIPushTab(v));
957:     PetscCall(PetscDualSpaceView(sum->sumspaces[i], v));
958:     PetscCall(PetscViewerASCIIPopTab(v));
959:   }
960:   PetscFunctionReturn(PETSC_SUCCESS);
961: }

963: static PetscErrorCode PetscDualSpaceView_Sum(PetscDualSpace sp, PetscViewer viewer)
964: {
965:   PetscBool iascii;

967:   PetscFunctionBegin;
968:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
969:   if (iascii) PetscCall(PetscDualSpaceSumView_Ascii(sp, viewer));
970:   PetscFunctionReturn(PETSC_SUCCESS);
971: }

973: static PetscErrorCode PetscDualSpaceDestroy_Sum(PetscDualSpace sp)
974: {
975:   PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;
976:   PetscInt            i, Ns = sum->numSumSpaces;

978:   PetscFunctionBegin;
979:   if (sum->symperms) {
980:     PetscInt **selfSyms = sum->symperms[0];

982:     if (selfSyms) {
983:       PetscInt i, **allocated = &selfSyms[-sum->selfSymOff];

985:       for (i = 0; i < sum->numSelfSym; i++) PetscCall(PetscFree(allocated[i]));
986:       PetscCall(PetscFree(allocated));
987:     }
988:     PetscCall(PetscFree(sum->symperms));
989:   }
990:   if (sum->symflips) {
991:     PetscScalar **selfSyms = sum->symflips[0];

993:     if (selfSyms) {
994:       PetscInt      i;
995:       PetscScalar **allocated = &selfSyms[-sum->selfSymOff];

997:       for (i = 0; i < sum->numSelfSym; i++) PetscCall(PetscFree(allocated[i]));
998:       PetscCall(PetscFree(allocated));
999:     }
1000:     PetscCall(PetscFree(sum->symflips));
1001:   }
1002:   for (i = 0; i < Ns; ++i) {
1003:     PetscCall(PetscDualSpaceDestroy(&sum->sumspaces[i]));
1004:     if (sum->all_rows) PetscCall(ISLocalToGlobalMappingDestroy(&sum->all_rows[i]));
1005:     if (sum->all_cols) PetscCall(ISLocalToGlobalMappingDestroy(&sum->all_cols[i]));
1006:     if (sum->int_rows) PetscCall(ISLocalToGlobalMappingDestroy(&sum->int_rows[i]));
1007:     if (sum->int_cols) PetscCall(ISLocalToGlobalMappingDestroy(&sum->int_cols[i]));
1008:   }
1009:   PetscCall(PetscFree(sum->sumspaces));
1010:   PetscCall(PetscFree(sum->all_rows));
1011:   PetscCall(PetscFree(sum->all_cols));
1012:   PetscCall(PetscFree(sum->int_rows));
1013:   PetscCall(PetscFree(sum->int_cols));
1014:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetSubspace_C", NULL));
1015:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetSubspace_C", NULL));
1016:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetNumSubspaces_C", NULL));
1017:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetNumSubspaces_C", NULL));
1018:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetConcatenate_C", NULL));
1019:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetConcatenate_C", NULL));
1020:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetInterleave_C", NULL));
1021:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetInterleave_C", NULL));
1022:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL));
1023:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL));
1024:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetMomentOrder_C", NULL));
1025:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetMomentOrder_C", NULL));
1026:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetNodeType_C", NULL));
1027:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetNodeType_C", NULL));
1028:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTensor_C", NULL));
1029:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTensor_C", NULL));
1030:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTrimmed_C", NULL));
1031:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTrimmed_C", NULL));
1032:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetUseMoments_C", NULL));
1033:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetUseMoments_C", NULL));
1034:   PetscCall(PetscFree(sum));
1035:   PetscFunctionReturn(PETSC_SUCCESS);
1036: }

1038: /*@
1039:   PetscDualSpaceSumSetInterleave - Set whether the basis functions and components of a uniform sum are interleaved

1041:   Logically collective

1043:   Input Parameters:
1044: + sp                    - a `PetscDualSpace` of type `PETSCDUALSPACESUM`
1045: . interleave_basis      - if `PETSC_TRUE`, the basis vectors of the subspaces are interleaved
1046: - interleave_components - if `PETSC_TRUE` and the space concatenates components (`PetscDualSpaceSumGetConcatenate()`),
1047:                           interleave the concatenated components

1049:   Level: developer

1051: .seealso: `PetscDualSpace`, `PETSCDUALSPACESUM`, `PETSCFEVECTOR`, `PetscDualSpaceSumGetInterleave()`
1052: @*/
1053: PetscErrorCode PetscDualSpaceSumSetInterleave(PetscDualSpace sp, PetscBool interleave_basis, PetscBool interleave_components)
1054: {
1055:   PetscFunctionBegin;
1057:   PetscTryMethod(sp, "PetscDualSpaceSumSetInterleave_C", (PetscDualSpace, PetscBool, PetscBool), (sp, interleave_basis, interleave_components));
1058:   PetscFunctionReturn(PETSC_SUCCESS);
1059: }

1061: static PetscErrorCode PetscDualSpaceSumSetInterleave_Sum(PetscDualSpace sp, PetscBool interleave_basis, PetscBool interleave_components)
1062: {
1063:   PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;

1065:   PetscFunctionBegin;
1066:   sum->interleave_basis      = interleave_basis;
1067:   sum->interleave_components = interleave_components;
1068:   PetscFunctionReturn(PETSC_SUCCESS);
1069: }

1071: /*@
1072:   PetscDualSpaceSumGetInterleave - Get whether the basis functions and components of a uniform sum are interleaved

1074:   Logically collective

1076:   Input Parameter:
1077: . sp - a `PetscDualSpace` of type `PETSCDUALSPACESUM`

1079:   Output Parameters:
1080: + interleave_basis      - if `PETSC_TRUE`, the basis vectors of the subspaces are interleaved
1081: - interleave_components - if `PETSC_TRUE` and the space concatenates components (`PetscDualSpaceSumGetConcatenate()`),
1082:                           interleave the concatenated components

1084:   Level: developer

1086: .seealso: `PetscDualSpace`, `PETSCDUALSPACESUM`, `PETSCFEVECTOR`, `PetscDualSpaceSumSetInterleave()`
1087: @*/
1088: PetscErrorCode PetscDualSpaceSumGetInterleave(PetscDualSpace sp, PetscBool *interleave_basis, PetscBool *interleave_components)
1089: {
1090:   PetscFunctionBegin;
1092:   if (interleave_basis) PetscAssertPointer(interleave_basis, 2);
1093:   if (interleave_components) PetscAssertPointer(interleave_components, 3);
1094:   PetscTryMethod(sp, "PetscDualSpaceSumGetInterleave_C", (PetscDualSpace, PetscBool *, PetscBool *), (sp, interleave_basis, interleave_components));
1095:   PetscFunctionReturn(PETSC_SUCCESS);
1096: }

1098: static PetscErrorCode PetscDualSpaceSumGetInterleave_Sum(PetscDualSpace sp, PetscBool *interleave_basis, PetscBool *interleave_components)
1099: {
1100:   PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data;

1102:   PetscFunctionBegin;
1103:   if (interleave_basis) *interleave_basis = sum->interleave_basis;
1104:   if (interleave_components) *interleave_components = sum->interleave_components;
1105:   PetscFunctionReturn(PETSC_SUCCESS);
1106: }

1108: #define PetscDualSpaceSumPassthrough(sp, func, ...) \
1109:   do { \
1110:     PetscDualSpace_Sum *sum = (PetscDualSpace_Sum *)sp->data; \
1111:     PetscBool           is_uniform; \
1112:     PetscCall(PetscDualSpaceSumIsUniform(sp, &is_uniform)); \
1113:     if (is_uniform && sum->numSumSpaces > 0) { \
1114:       PetscDualSpace subsp; \
1115:       PetscCall(PetscDualSpaceSumGetSubspace(sp, 0, &subsp)); \
1116:       PetscCall(func(subsp, __VA_ARGS__)); \
1117:     } \
1118:   } while (0)

1120: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Sum(PetscDualSpace sp, PetscBool *value)
1121: {
1122:   PetscFunctionBegin;
1123:   PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetContinuity, value);
1124:   PetscFunctionReturn(PETSC_SUCCESS);
1125: }

1127: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Sum(PetscDualSpace sp, PetscBool value)
1128: {
1129:   PetscFunctionBegin;
1130:   PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetContinuity, value);
1131:   PetscFunctionReturn(PETSC_SUCCESS);
1132: }

1134: static PetscErrorCode PetscDualSpaceLagrangeGetTensor_Sum(PetscDualSpace sp, PetscBool *value)
1135: {
1136:   PetscFunctionBegin;
1137:   PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetTensor, value);
1138:   PetscFunctionReturn(PETSC_SUCCESS);
1139: }

1141: static PetscErrorCode PetscDualSpaceLagrangeSetTensor_Sum(PetscDualSpace sp, PetscBool value)
1142: {
1143:   PetscFunctionBegin;
1144:   PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetTensor, value);
1145:   PetscFunctionReturn(PETSC_SUCCESS);
1146: }

1148: static PetscErrorCode PetscDualSpaceLagrangeGetTrimmed_Sum(PetscDualSpace sp, PetscBool *value)
1149: {
1150:   PetscFunctionBegin;
1151:   PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetTrimmed, value);
1152:   PetscFunctionReturn(PETSC_SUCCESS);
1153: }

1155: static PetscErrorCode PetscDualSpaceLagrangeSetTrimmed_Sum(PetscDualSpace sp, PetscBool value)
1156: {
1157:   PetscFunctionBegin;
1158:   PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetTrimmed, value);
1159:   PetscFunctionReturn(PETSC_SUCCESS);
1160: }

1162: static PetscErrorCode PetscDualSpaceLagrangeGetUseMoments_Sum(PetscDualSpace sp, PetscBool *value)
1163: {
1164:   PetscFunctionBegin;
1165:   PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetUseMoments, value);
1166:   PetscFunctionReturn(PETSC_SUCCESS);
1167: }

1169: static PetscErrorCode PetscDualSpaceLagrangeSetUseMoments_Sum(PetscDualSpace sp, PetscBool value)
1170: {
1171:   PetscFunctionBegin;
1172:   PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetUseMoments, value);
1173:   PetscFunctionReturn(PETSC_SUCCESS);
1174: }

1176: static PetscErrorCode PetscDualSpaceLagrangeGetMomentOrder_Sum(PetscDualSpace sp, PetscInt *value)
1177: {
1178:   PetscFunctionBegin;
1179:   PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetMomentOrder, value);
1180:   PetscFunctionReturn(PETSC_SUCCESS);
1181: }

1183: static PetscErrorCode PetscDualSpaceLagrangeSetMomentOrder_Sum(PetscDualSpace sp, PetscInt value)
1184: {
1185:   PetscFunctionBegin;
1186:   PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetMomentOrder, value);
1187:   PetscFunctionReturn(PETSC_SUCCESS);
1188: }

1190: static PetscErrorCode PetscDualSpaceLagrangeGetNodeType_Sum(PetscDualSpace sp, PetscDTNodeType *node_type, PetscBool *include_endpoints, PetscReal *exponent)
1191: {
1192:   PetscFunctionBegin;
1193:   PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeGetNodeType, node_type, include_endpoints, exponent);
1194:   PetscFunctionReturn(PETSC_SUCCESS);
1195: }

1197: static PetscErrorCode PetscDualSpaceLagrangeSetNodeType_Sum(PetscDualSpace sp, PetscDTNodeType node_type, PetscBool include_endpoints, PetscReal exponent)
1198: {
1199:   PetscFunctionBegin;
1200:   PetscDualSpaceSumPassthrough(sp, PetscDualSpaceLagrangeSetNodeType, node_type, include_endpoints, exponent);
1201:   PetscFunctionReturn(PETSC_SUCCESS);
1202: }

1204: static PetscErrorCode PetscDualSpaceInitialize_Sum(PetscDualSpace sp)
1205: {
1206:   PetscFunctionBegin;
1207:   sp->ops->destroy              = PetscDualSpaceDestroy_Sum;
1208:   sp->ops->view                 = PetscDualSpaceView_Sum;
1209:   sp->ops->setfromoptions       = NULL;
1210:   sp->ops->duplicate            = PetscDualSpaceDuplicate_Sum;
1211:   sp->ops->setup                = PetscDualSpaceSetUp_Sum;
1212:   sp->ops->createheightsubspace = NULL;
1213:   sp->ops->createpointsubspace  = NULL;
1214:   sp->ops->getsymmetries        = PetscDualSpaceGetSymmetries_Sum;
1215:   sp->ops->apply                = PetscDualSpaceApplyDefault;
1216:   sp->ops->applyall             = PetscDualSpaceApplyAllDefault;
1217:   sp->ops->applyint             = PetscDualSpaceApplyInteriorDefault;
1218:   sp->ops->createalldata        = PetscDualSpaceCreateAllDataDefault;
1219:   sp->ops->createintdata        = PetscDualSpaceCreateInteriorDataDefault;

1221:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetNumSubspaces_C", PetscDualSpaceSumGetNumSubspaces_Sum));
1222:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetNumSubspaces_C", PetscDualSpaceSumSetNumSubspaces_Sum));
1223:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetSubspace_C", PetscDualSpaceSumGetSubspace_Sum));
1224:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetSubspace_C", PetscDualSpaceSumSetSubspace_Sum));
1225:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetConcatenate_C", PetscDualSpaceSumGetConcatenate_Sum));
1226:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetConcatenate_C", PetscDualSpaceSumSetConcatenate_Sum));
1227:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumGetInterleave_C", PetscDualSpaceSumGetInterleave_Sum));
1228:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceSumSetInterleave_C", PetscDualSpaceSumSetInterleave_Sum));
1229:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Sum));
1230:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Sum));
1231:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetMomentOrder_C", PetscDualSpaceLagrangeGetMomentOrder_Sum));
1232:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetMomentOrder_C", PetscDualSpaceLagrangeSetMomentOrder_Sum));
1233:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetNodeType_C", PetscDualSpaceLagrangeGetNodeType_Sum));
1234:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetNodeType_C", PetscDualSpaceLagrangeSetNodeType_Sum));
1235:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTensor_C", PetscDualSpaceLagrangeGetTensor_Sum));
1236:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTensor_C", PetscDualSpaceLagrangeSetTensor_Sum));
1237:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetTrimmed_C", PetscDualSpaceLagrangeGetTrimmed_Sum));
1238:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetTrimmed_C", PetscDualSpaceLagrangeSetTrimmed_Sum));
1239:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeGetUseMoments_C", PetscDualSpaceLagrangeGetUseMoments_Sum));
1240:   PetscCall(PetscObjectComposeFunction((PetscObject)sp, "PetscDualSpaceLagrangeSetUseMoments_C", PetscDualSpaceLagrangeSetUseMoments_Sum));
1241:   PetscFunctionReturn(PETSC_SUCCESS);
1242: }

1244: /*MC
1245:   PETSCDUALSPACESUM = "sum" - A `PetscDualSpace` object that encapsulates a sum of subspaces.

1247:   Level: intermediate

1249:   Note:
1250:   That sum can either be direct or a concatenation. For example if A and B are spaces each with 2 components,
1251:   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
1252:   same reference element.

1254: .seealso: `PetscDualSpace`, `PetscDualSpaceType`, `PetscDualSpaceCreate()`, `PetscDualSpaceSetType()`, `PetscDualSpaceSumGetNumSubspaces()`, `PetscDualSpaceSumSetNumSubspaces()`,
1255:           `PetscDualSpaceSumGetConcatenate()`, `PetscDualSpaceSumSetConcatenate()`, `PetscDualSpaceSumSetInterleave()`, `PetscDualSpaceSumGetInterleave()`
1256: M*/
1257: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Sum(PetscDualSpace sp)
1258: {
1259:   PetscDualSpace_Sum *sum;

1261:   PetscFunctionBegin;
1263:   PetscCall(PetscNew(&sum));
1264:   sum->numSumSpaces = PETSC_DEFAULT;
1265:   sp->data          = sum;
1266:   sp->k             = PETSC_FORM_DEGREE_UNDEFINED;
1267:   PetscCall(PetscDualSpaceInitialize_Sum(sp));
1268:   PetscFunctionReturn(PETSC_SUCCESS);
1269: }

1271: /*@
1272:   PetscDualSpaceCreateSum - Create a finite element dual basis that is the sum of other dual bases

1274:   Collective

1276:   Input Parameters:
1277: + numSubspaces - the number of spaces that will be added together
1278: . subspaces    - an array of length `numSubspaces` of spaces
1279: - concatenate  - if `PETSC_FALSE`, the sum-space has the same components as the individual dual spaces (`PetscDualSpaceGetNumComponents()`); if `PETSC_TRUE`, the individual components are concatenated to create a dual space with more components

1281:   Output Parameter:
1282: . sumSpace - a `PetscDualSpace` of type `PETSCDUALSPACESUM`

1284:   Level: advanced

1286: .seealso: `PetscDualSpace`, `PETSCDUALSPACESUM`, `PETSCSPACESUM`
1287: @*/
1288: PetscErrorCode PetscDualSpaceCreateSum(PetscInt numSubspaces, const PetscDualSpace subspaces[], PetscBool concatenate, PetscDualSpace *sumSpace)
1289: {
1290:   PetscInt i, Nc = 0;

1292:   PetscFunctionBegin;
1293:   PetscCall(PetscDualSpaceCreate(PetscObjectComm((PetscObject)subspaces[0]), sumSpace));
1294:   PetscCall(PetscDualSpaceSetType(*sumSpace, PETSCDUALSPACESUM));
1295:   PetscCall(PetscDualSpaceSumSetNumSubspaces(*sumSpace, numSubspaces));
1296:   PetscCall(PetscDualSpaceSumSetConcatenate(*sumSpace, concatenate));
1297:   for (i = 0; i < numSubspaces; ++i) {
1298:     PetscInt sNc;

1300:     PetscCall(PetscDualSpaceSumSetSubspace(*sumSpace, i, subspaces[i]));
1301:     PetscCall(PetscDualSpaceGetNumComponents(subspaces[i], &sNc));
1302:     if (concatenate) Nc += sNc;
1303:     else Nc = sNc;

1305:     if (i == 0) {
1306:       DM dm;

1308:       PetscCall(PetscDualSpaceGetDM(subspaces[i], &dm));
1309:       PetscCall(PetscDualSpaceSetDM(*sumSpace, dm));
1310:     }
1311:   }
1312:   PetscCall(PetscDualSpaceSetNumComponents(*sumSpace, Nc));
1313:   PetscFunctionReturn(PETSC_SUCCESS);
1314: }