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 isascii;

    967:   PetscFunctionBegin;
    968:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    969:   if (isascii) 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, PeOp PetscBool *interleave_basis, PeOp 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: }