Actual source code: dmi.c

      1: #include <petsc/private/dmimpl.h>
      2: #include <petscds.h>

      4: // Greatest common divisor of two nonnegative integers
      5: PetscInt PetscGCD(PetscInt a, PetscInt b)
      6: {
      7:   while (b != 0) {
      8:     PetscInt tmp = b;
      9:     b            = a % b;
     10:     a            = tmp;
     11:   }
     12:   return a;
     13: }

     15: PetscErrorCode DMCreateGlobalVector_Section_Private(DM dm, Vec *vec)
     16: {
     17:   PetscSection gSection;
     18:   PetscInt     localSize, bs, blockSize = -1, pStart, pEnd, p;
     19:   PetscInt     in[2], out[2];

     21:   PetscFunctionBegin;
     22:   PetscCall(DMGetGlobalSection(dm, &gSection));
     23:   PetscCall(PetscSectionGetChart(gSection, &pStart, &pEnd));
     24:   for (p = pStart; p < pEnd; ++p) {
     25:     PetscInt dof, cdof;

     27:     PetscCall(PetscSectionGetDof(gSection, p, &dof));
     28:     PetscCall(PetscSectionGetConstraintDof(gSection, p, &cdof));

     30:     if (dof - cdof > 0) {
     31:       if (blockSize < 0) {
     32:         /* set blockSize */
     33:         blockSize = dof - cdof;
     34:       } else {
     35:         blockSize = PetscGCD(dof - cdof, blockSize);
     36:       }
     37:     }
     38:   }

     40:   // You cannot negate PETSC_INT_MIN
     41:   in[0] = blockSize < 0 ? -PETSC_INT_MAX : -blockSize;
     42:   in[1] = blockSize;
     43:   PetscCallMPI(MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm)));
     44:   /* -out[0] = min(blockSize), out[1] = max(blockSize) */
     45:   if (-out[0] == out[1]) {
     46:     bs = out[1];
     47:   } else bs = 1;

     49:   if (bs < 0) { /* Everyone was empty */
     50:     blockSize = 1;
     51:     bs        = 1;
     52:   }

     54:   PetscCall(PetscSectionGetConstrainedStorageSize(gSection, &localSize));
     55:   PetscCheck(localSize % blockSize == 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mismatch between blocksize %" PetscInt_FMT " and local storage size %" PetscInt_FMT, blockSize, localSize);
     56:   PetscCall(VecCreate(PetscObjectComm((PetscObject)dm), vec));
     57:   PetscCall(VecSetSizes(*vec, localSize, PETSC_DETERMINE));
     58:   PetscCall(VecSetBlockSize(*vec, bs));
     59:   PetscCall(VecSetType(*vec, dm->vectype));
     60:   PetscCall(VecSetDM(*vec, dm));
     61:   /* PetscCall(VecSetLocalToGlobalMapping(*vec, dm->ltogmap)); */
     62:   PetscFunctionReturn(PETSC_SUCCESS);
     63: }

     65: PetscErrorCode DMCreateLocalVector_Section_Private(DM dm, Vec *vec)
     66: {
     67:   PetscSection section;
     68:   PetscInt     localSize, blockSize = -1, pStart, pEnd, p;

     70:   PetscFunctionBegin;
     71:   PetscCall(DMGetLocalSection(dm, &section));
     72:   PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
     73:   for (p = pStart; p < pEnd; ++p) {
     74:     PetscInt dof;

     76:     PetscCall(PetscSectionGetDof(section, p, &dof));
     77:     if ((blockSize < 0) && (dof > 0)) blockSize = dof;
     78:     if (dof > 0) blockSize = PetscGCD(dof, blockSize);
     79:   }
     80:   PetscCall(PetscSectionGetStorageSize(section, &localSize));
     81:   PetscCall(VecCreate(PETSC_COMM_SELF, vec));
     82:   PetscCall(VecSetSizes(*vec, localSize, localSize));
     83:   PetscCall(VecSetBlockSize(*vec, blockSize));
     84:   PetscCall(VecSetType(*vec, dm->vectype));
     85:   PetscCall(VecSetDM(*vec, dm));
     86:   PetscFunctionReturn(PETSC_SUCCESS);
     87: }

     89: static PetscErrorCode PetscSectionSelectFields_Private(PetscSection s, PetscSection gs, PetscInt numFields, const PetscInt fields[], const PetscInt numComps[], const PetscInt comps[], IS *is)
     90: {
     91:   PetscInt *subIndices;
     92:   PetscInt  bs = 0, bsLocal[2], bsMinMax[2];
     93:   PetscInt  pStart, pEnd, Nc, subSize = 0, subOff = 0;

     95:   PetscFunctionBegin;
     96:   PetscCall(PetscSectionGetChart(gs, &pStart, &pEnd));
     97:   if (numComps) {
     98:     for (PetscInt f = 0, off = 0; f < numFields; ++f) {
     99:       PetscInt Nc;

    101:       if (numComps[f] < 0) continue;
    102:       PetscCall(PetscSectionGetFieldComponents(s, f, &Nc));
    103:       PetscCheck(numComps[f] <= Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT ": Number of selected components %" PetscInt_FMT " > %" PetscInt_FMT " number of field components", f, numComps[f], Nc);
    104:       for (PetscInt c = 0; c < numComps[f]; ++c, ++off) PetscCheck(comps[off] < Nc, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field %" PetscInt_FMT ": component %" PetscInt_FMT " not in [0, %" PetscInt_FMT ")", f, comps[off], Nc);
    105:       bs += numComps[f];
    106:     }
    107:   } else {
    108:     for (PetscInt f = 0; f < numFields; ++f) {
    109:       PetscInt Nc;

    111:       PetscCall(PetscSectionGetFieldComponents(s, fields[f], &Nc));
    112:       bs += Nc;
    113:     }
    114:   }
    115:   for (PetscInt p = pStart; p < pEnd; ++p) {
    116:     PetscInt gdof, pSubSize = 0;

    118:     PetscCall(PetscSectionGetDof(gs, p, &gdof));
    119:     if (gdof > 0) {
    120:       PetscInt off = 0;

    122:       for (PetscInt f = 0; f < numFields; ++f) {
    123:         PetscInt fdof, fcdof, sfdof, sfcdof = 0;

    125:         PetscCall(PetscSectionGetFieldComponents(s, f, &Nc));
    126:         PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
    127:         PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &fcdof));
    128:         if (numComps && numComps[f] >= 0) {
    129:           const PetscInt *ind;

    131:           // Assume sets of dofs on points are of size Nc
    132:           PetscCheck(!(fdof % Nc), PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of components %" PetscInt_FMT " should evenly divide the dofs %" PetscInt_FMT " on point %" PetscInt_FMT, Nc, fdof, p);
    133:           sfdof = (fdof / Nc) * numComps[f];
    134:           PetscCall(PetscSectionGetFieldConstraintIndices(s, p, fields[f], &ind));
    135:           for (PetscInt i = 0; i < (fdof / Nc); ++i) {
    136:             for (PetscInt c = 0, fcc = 0; c < Nc; ++c) {
    137:               if (c == comps[off + fcc]) {
    138:                 ++fcc;
    139:                 ++sfcdof;
    140:               }
    141:             }
    142:           }
    143:           pSubSize += sfdof - sfcdof;
    144:           off += numComps[f];
    145:         } else {
    146:           pSubSize += fdof - fcdof;
    147:         }
    148:       }
    149:       subSize += pSubSize;
    150:       if (pSubSize && bs != pSubSize) {
    151:         // Layout does not admit a pointwise block size
    152:         bs = 1;
    153:       }
    154:     }
    155:   }
    156:   // Must have same blocksize on all procs (some might have no points)
    157:   bsLocal[0] = bs < 0 ? PETSC_INT_MAX : bs;
    158:   bsLocal[1] = bs;
    159:   PetscCall(PetscGlobalMinMaxInt(PetscObjectComm((PetscObject)gs), bsLocal, bsMinMax));
    160:   if (bsMinMax[0] != bsMinMax[1]) {
    161:     bs = 1;
    162:   } else {
    163:     bs = bsMinMax[0];
    164:   }
    165:   PetscCall(PetscMalloc1(subSize, &subIndices));
    166:   for (PetscInt p = pStart; p < pEnd; ++p) {
    167:     PetscInt gdof, goff;

    169:     PetscCall(PetscSectionGetDof(gs, p, &gdof));
    170:     if (gdof > 0) {
    171:       PetscInt off = 0;

    173:       PetscCall(PetscSectionGetOffset(gs, p, &goff));
    174:       for (PetscInt f = 0; f < numFields; ++f) {
    175:         PetscInt fdof, fcdof, poff = 0;

    177:         /* Can get rid of this loop by storing field information in the global section */
    178:         for (PetscInt f2 = 0; f2 < fields[f]; ++f2) {
    179:           PetscCall(PetscSectionGetFieldDof(s, p, f2, &fdof));
    180:           PetscCall(PetscSectionGetFieldConstraintDof(s, p, f2, &fcdof));
    181:           poff += fdof - fcdof;
    182:         }
    183:         PetscCall(PetscSectionGetFieldDof(s, p, fields[f], &fdof));
    184:         PetscCall(PetscSectionGetFieldConstraintDof(s, p, fields[f], &fcdof));

    186:         if (numComps && numComps[f] >= 0) {
    187:           const PetscInt *ind;

    189:           // Assume sets of dofs on points are of size Nc
    190:           PetscCall(PetscSectionGetFieldConstraintIndices(s, p, fields[f], &ind));
    191:           for (PetscInt i = 0, fcoff = 0, pfoff = 0; i < (fdof / Nc); ++i) {
    192:             for (PetscInt c = 0, fcc = 0; c < Nc; ++c) {
    193:               const PetscInt k = i * Nc + c;

    195:               if (ind[fcoff] == k) {
    196:                 ++fcoff;
    197:                 continue;
    198:               }
    199:               if (c == comps[off + fcc]) {
    200:                 ++fcc;
    201:                 subIndices[subOff++] = goff + poff + pfoff;
    202:               }
    203:               ++pfoff;
    204:             }
    205:           }
    206:           off += numComps[f];
    207:         } else {
    208:           for (PetscInt fc = 0; fc < fdof - fcdof; ++fc, ++subOff) subIndices[subOff] = goff + poff + fc;
    209:         }
    210:       }
    211:     }
    212:   }
    213:   PetscCheck(subSize == subOff, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "The offset array size %" PetscInt_FMT " != %" PetscInt_FMT " the number of indices", subSize, subOff);
    214:   PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)gs), subSize, subIndices, PETSC_OWN_POINTER, is));
    215:   if (bs > 1) {
    216:     // We need to check that the block size does not come from non-contiguous fields
    217:     PetscInt set = 1, rset = 1;
    218:     for (PetscInt i = 0; i < subSize; i += bs) {
    219:       for (PetscInt j = 0; j < bs; ++j) {
    220:         if (subIndices[i + j] != subIndices[i] + j) {
    221:           set = 0;
    222:           break;
    223:         }
    224:       }
    225:     }
    226:     PetscCallMPI(MPIU_Allreduce(&set, &rset, 1, MPIU_INT, MPI_PROD, PetscObjectComm((PetscObject)gs)));
    227:     if (rset) PetscCall(ISSetBlockSize(*is, bs));
    228:   }
    229:   PetscFunctionReturn(PETSC_SUCCESS);
    230: }

    232: static PetscErrorCode DMSelectFields_Private(DM dm, PetscSection section, PetscInt numFields, const PetscInt fields[], const PetscInt numComps[], const PetscInt comps[], IS *is, DM *subdm)
    233: {
    234:   PetscSection subsection;
    235:   PetscBool    haveNull = PETSC_FALSE;
    236:   PetscInt     nf = 0, of = 0;

    238:   PetscFunctionBegin;
    239:   if (numComps) {
    240:     const PetscInt field = fields[0];

    242:     PetscCheck(numFields == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "We only support a single field for component selection right now");
    243:     PetscCall(PetscSectionCreateComponentSubsection(section, numComps[field], comps, &subsection));
    244:     PetscCall(DMSetLocalSection(*subdm, subsection));
    245:     PetscCall(PetscSectionDestroy(&subsection));
    246:     (*subdm)->nullspaceConstructors[field] = dm->nullspaceConstructors[field];
    247:     if (dm->probs) {
    248:       PetscFV  fv, fvNew;
    249:       PetscInt fnum[1] = {field};

    251:       PetscCall(DMSetNumFields(*subdm, 1));
    252:       PetscCall(DMGetField(dm, field, NULL, (PetscObject *)&fv));
    253:       PetscCall(PetscFVClone(fv, &fvNew));
    254:       PetscCall(PetscFVSetNumComponents(fvNew, numComps[0]));
    255:       PetscCall(DMSetField(*subdm, 0, NULL, (PetscObject)fvNew));
    256:       PetscCall(PetscFVDestroy(&fvNew));
    257:       PetscCall(DMCreateDS(*subdm));
    258:       if (numComps[0] == 1 && is) {
    259:         PetscObject disc, space, pmat;

    261:         PetscCall(DMGetField(*subdm, field, NULL, &disc));
    262:         PetscCall(PetscObjectQuery(disc, "nullspace", &space));
    263:         if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", space));
    264:         PetscCall(PetscObjectQuery(disc, "nearnullspace", &space));
    265:         if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nearnullspace", space));
    266:         PetscCall(PetscObjectQuery(disc, "pmat", &pmat));
    267:         if (pmat) PetscCall(PetscObjectCompose((PetscObject)*is, "pmat", pmat));
    268:       }
    269:       PetscCall(PetscDSCopyConstants(dm->probs[field].ds, (*subdm)->probs[0].ds));
    270:       PetscCall(PetscDSCopyBoundary(dm->probs[field].ds, 1, fnum, (*subdm)->probs[0].ds));
    271:       PetscCall(PetscDSSelectEquations(dm->probs[field].ds, 1, fnum, (*subdm)->probs[0].ds));
    272:     }
    273:     if ((*subdm)->nullspaceConstructors[0] && is) {
    274:       MatNullSpace nullSpace;

    276:       PetscCall((*(*subdm)->nullspaceConstructors[0])(*subdm, 0, 0, &nullSpace));
    277:       PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", (PetscObject)nullSpace));
    278:       PetscCall(MatNullSpaceDestroy(&nullSpace));
    279:     }
    280:     if (dm->coarseMesh) PetscCall(DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh));
    281:     PetscFunctionReturn(PETSC_SUCCESS);
    282:   }
    283:   PetscCall(PetscSectionCreateSubsection(section, numFields, fields, &subsection));
    284:   PetscCall(DMSetLocalSection(*subdm, subsection));
    285:   PetscCall(PetscSectionDestroy(&subsection));
    286:   for (PetscInt f = 0; f < numFields; ++f) {
    287:     (*subdm)->nullspaceConstructors[f] = dm->nullspaceConstructors[fields[f]];
    288:     if ((*subdm)->nullspaceConstructors[f]) {
    289:       haveNull = PETSC_TRUE;
    290:       nf       = f;
    291:       of       = fields[f];
    292:     }
    293:   }
    294:   if (dm->probs) {
    295:     PetscCall(DMSetNumFields(*subdm, numFields));
    296:     for (PetscInt f = 0; f < numFields; ++f) {
    297:       PetscObject disc;

    299:       PetscCall(DMGetField(dm, fields[f], NULL, &disc));
    300:       PetscCall(DMSetField(*subdm, f, NULL, disc));
    301:     }
    302:     // TODO: if only FV, then cut down the components
    303:     PetscCall(DMCreateDS(*subdm));
    304:     if (numFields == 1 && is) {
    305:       PetscObject disc, space, pmat;

    307:       PetscCall(DMGetField(*subdm, 0, NULL, &disc));
    308:       PetscCall(PetscObjectQuery(disc, "nullspace", &space));
    309:       if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", space));
    310:       PetscCall(PetscObjectQuery(disc, "nearnullspace", &space));
    311:       if (space) PetscCall(PetscObjectCompose((PetscObject)*is, "nearnullspace", space));
    312:       PetscCall(PetscObjectQuery(disc, "pmat", &pmat));
    313:       if (pmat) PetscCall(PetscObjectCompose((PetscObject)*is, "pmat", pmat));
    314:     }
    315:     // Check if DSes record their DM fields
    316:     if (dm->probs[0].fields) {
    317:       PetscInt d, e;

    319:       for (d = 0, e = 0; d < dm->Nds && e < (*subdm)->Nds; ++d) {
    320:         const PetscInt  Nf = dm->probs[d].ds->Nf;
    321:         const PetscInt *fld;
    322:         PetscInt        f, g;

    324:         PetscCall(ISGetIndices(dm->probs[d].fields, &fld));
    325:         for (f = 0; f < Nf; ++f) {
    326:           for (g = 0; g < numFields; ++g)
    327:             if (fld[f] == fields[g]) break;
    328:           if (g < numFields) break;
    329:         }
    330:         PetscCall(ISRestoreIndices(dm->probs[d].fields, &fld));
    331:         if (f == Nf) continue;
    332:         PetscCall(PetscDSCopyConstants(dm->probs[d].ds, (*subdm)->probs[e].ds));
    333:         PetscCall(PetscDSCopyBoundary(dm->probs[d].ds, numFields, fields, (*subdm)->probs[e].ds));
    334:         // Translate DM fields to DS fields
    335:         {
    336:           IS              infields, dsfields;
    337:           const PetscInt *fld, *ofld;
    338:           PetscInt       *fidx;
    339:           PetscInt        onf, nf;

    341:           PetscCall(ISCreateGeneral(PETSC_COMM_SELF, numFields, fields, PETSC_USE_POINTER, &infields));
    342:           PetscCall(ISIntersect(infields, dm->probs[d].fields, &dsfields));
    343:           PetscCall(ISDestroy(&infields));
    344:           PetscCall(ISGetLocalSize(dsfields, &nf));
    345:           PetscCheck(nf, PETSC_COMM_SELF, PETSC_ERR_PLIB, "DS cannot be supported on 0 fields");
    346:           PetscCall(ISGetIndices(dsfields, &fld));
    347:           PetscCall(ISGetLocalSize(dm->probs[d].fields, &onf));
    348:           PetscCall(ISGetIndices(dm->probs[d].fields, &ofld));
    349:           PetscCall(PetscMalloc1(nf, &fidx));
    350:           for (PetscInt f = 0, g = 0; f < onf && g < nf; ++f) {
    351:             if (ofld[f] == fld[g]) fidx[g++] = f;
    352:           }
    353:           PetscCall(ISRestoreIndices(dm->probs[d].fields, &ofld));
    354:           PetscCall(ISRestoreIndices(dsfields, &fld));
    355:           PetscCall(ISDestroy(&dsfields));
    356:           PetscCall(PetscDSSelectDiscretizations(dm->probs[0].ds, nf, fidx, PETSC_DETERMINE, PETSC_DETERMINE, (*subdm)->probs[0].ds));
    357:           PetscCall(PetscDSSelectEquations(dm->probs[0].ds, nf, fidx, (*subdm)->probs[0].ds));
    358:           PetscCall(PetscFree(fidx));
    359:         }
    360:         ++e;
    361:       }
    362:     } else {
    363:       PetscCall(PetscDSCopyConstants(dm->probs[0].ds, (*subdm)->probs[0].ds));
    364:       PetscCall(PetscDSCopyBoundary(dm->probs[0].ds, PETSC_DETERMINE, NULL, (*subdm)->probs[0].ds));
    365:       PetscCall(PetscDSSelectDiscretizations(dm->probs[0].ds, numFields, fields, PETSC_DETERMINE, PETSC_DETERMINE, (*subdm)->probs[0].ds));
    366:       PetscCall(PetscDSSelectEquations(dm->probs[0].ds, numFields, fields, (*subdm)->probs[0].ds));
    367:     }
    368:   }
    369:   if (haveNull && is) {
    370:     MatNullSpace nullSpace;

    372:     PetscCall((*(*subdm)->nullspaceConstructors[nf])(*subdm, of, nf, &nullSpace));
    373:     PetscCall(PetscObjectCompose((PetscObject)*is, "nullspace", (PetscObject)nullSpace));
    374:     PetscCall(MatNullSpaceDestroy(&nullSpace));
    375:   }
    376:   if (dm->coarseMesh) PetscCall(DMCreateSubDM(dm->coarseMesh, numFields, fields, NULL, &(*subdm)->coarseMesh));
    377:   PetscFunctionReturn(PETSC_SUCCESS);
    378: }

    380: /*@
    381:   DMCreateSectionSubDM - Returns an `IS` and `subDM` containing a `PetscSection` that encapsulates a subproblem defined by a subset of the fields in a `PetscSection` in the `DM`.

    383:   Not Collective

    385:   Input Parameters:
    386: + dm        - The `DM` object
    387: . numFields - The number of fields to incorporate into `subdm`
    388: . fields    - The field numbers of the selected fields
    389: . numComps  - The number of components from each field to incorporate into `subdm`, or PETSC_DECIDE for all components
    390: - comps     - The component numbers of the selected fields (omitted for PTESC_DECIDE fields)

    392:   Output Parameters:
    393: + is    - The global indices for the subproblem or `NULL`
    394: - subdm - The `DM` for the subproblem, which must already have be cloned from `dm` or `NULL`

    396:   Level: intermediate

    398:   Notes:
    399:   If `is` and `subdm` are both `NULL` this does nothing

    401: .seealso: `DMCreateSubDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()`
    402: @*/
    403: PetscErrorCode DMCreateSectionSubDM(DM dm, PetscInt numFields, const PetscInt fields[], const PetscInt numComps[], const PetscInt comps[], IS *is, DM *subdm)
    404: {
    405:   PetscSection section, sectionGlobal;
    406:   PetscInt     Nf;

    408:   PetscFunctionBegin;
    409:   if (!numFields) PetscFunctionReturn(PETSC_SUCCESS);
    410:   PetscCall(DMGetLocalSection(dm, &section));
    411:   PetscCall(DMGetGlobalSection(dm, &sectionGlobal));
    412:   PetscCheck(section, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields");
    413:   PetscCheck(sectionGlobal, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields");
    414:   PetscCall(PetscSectionGetNumFields(section, &Nf));
    415:   PetscCheck(numFields <= Nf, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Number of requested fields %" PetscInt_FMT " greater than number of DM fields %" PetscInt_FMT, numFields, Nf);

    417:   if (is) PetscCall(PetscSectionSelectFields_Private(section, sectionGlobal, numFields, fields, numComps, comps, is));
    418:   if (subdm) PetscCall(DMSelectFields_Private(dm, section, numFields, fields, numComps, comps, is, subdm));
    419:   PetscFunctionReturn(PETSC_SUCCESS);
    420: }

    422: /*@C
    423:   DMCreateSectionSuperDM - Returns an arrays of `IS` and a `DM` containing a `PetscSection` that encapsulates a superproblem defined by the array of `DM` and their `PetscSection`

    425:   Not Collective

    427:   Input Parameters:
    428: + dms - The `DM` objects, the must all have the same topology; for example obtained with `DMClone()`
    429: - len - The number of `DM` in `dms`

    431:   Output Parameters:
    432: + is      - The global indices for the subproblem, or `NULL`
    433: - superdm - The `DM` for the superproblem, which must already have be cloned and contain the same topology as the `dms`

    435:   Level: intermediate

    437: .seealso: `DMCreateSuperDM()`, `DMGetLocalSection()`, `DMPlexSetMigrationSF()`, `DMView()`
    438: @*/
    439: PetscErrorCode DMCreateSectionSuperDM(DM dms[], PetscInt len, IS *is[], DM *superdm)
    440: {
    441:   MPI_Comm     comm;
    442:   PetscSection supersection, *sections, *sectionGlobals;
    443:   PetscInt    *Nfs, Nf = 0, f, supf, oldf = -1, nullf = -1, i;
    444:   PetscBool    haveNull = PETSC_FALSE;

    446:   PetscFunctionBegin;
    447:   PetscCall(PetscObjectGetComm((PetscObject)dms[0], &comm));
    448:   /* Pull out local and global sections */
    449:   PetscCall(PetscMalloc3(len, &Nfs, len, &sections, len, &sectionGlobals));
    450:   for (i = 0; i < len; ++i) {
    451:     PetscCall(DMGetLocalSection(dms[i], &sections[i]));
    452:     PetscCall(DMGetGlobalSection(dms[i], &sectionGlobals[i]));
    453:     PetscCheck(sections[i], comm, PETSC_ERR_ARG_WRONG, "Must set default section for DM before splitting fields");
    454:     PetscCheck(sectionGlobals[i], comm, PETSC_ERR_ARG_WRONG, "Must set default global section for DM before splitting fields");
    455:     PetscCall(PetscSectionGetNumFields(sections[i], &Nfs[i]));
    456:     Nf += Nfs[i];
    457:   }
    458:   /* Create the supersection */
    459:   PetscCall(PetscSectionCreateSupersection(sections, len, &supersection));
    460:   PetscCall(DMSetLocalSection(*superdm, supersection));
    461:   /* Create ISes */
    462:   if (is) {
    463:     PetscSection supersectionGlobal;
    464:     PetscInt     bs = -1, startf = 0;

    466:     PetscCall(PetscMalloc1(len, is));
    467:     PetscCall(DMGetGlobalSection(*superdm, &supersectionGlobal));
    468:     for (i = 0; i < len; startf += Nfs[i], ++i) {
    469:       PetscInt *subIndices;
    470:       PetscInt  subSize, subOff, pStart, pEnd, p, start, end, dummy;

    472:       PetscCall(PetscSectionGetChart(sectionGlobals[i], &pStart, &pEnd));
    473:       PetscCall(PetscSectionGetConstrainedStorageSize(sectionGlobals[i], &subSize));
    474:       PetscCall(PetscMalloc1(subSize, &subIndices));
    475:       for (p = pStart, subOff = 0; p < pEnd; ++p) {
    476:         PetscInt gdof, gcdof, gtdof, d;

    478:         PetscCall(PetscSectionGetDof(sectionGlobals[i], p, &gdof));
    479:         PetscCall(PetscSectionGetConstraintDof(sections[i], p, &gcdof));
    480:         gtdof = gdof - gcdof;
    481:         if (gdof > 0 && gtdof) {
    482:           if (bs < 0) {
    483:             bs = gtdof;
    484:           } else if (bs != gtdof) {
    485:             bs = 1;
    486:           }
    487:           PetscCall(DMGetGlobalFieldOffset_Private(*superdm, p, startf, &start, &dummy));
    488:           PetscCall(DMGetGlobalFieldOffset_Private(*superdm, p, startf + Nfs[i] - 1, &dummy, &end));
    489:           PetscCheck(end - start == gtdof, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid number of global dofs %" PetscInt_FMT " != %" PetscInt_FMT " for dm %" PetscInt_FMT " on point %" PetscInt_FMT, end - start, gtdof, i, p);
    490:           for (d = start; d < end; ++d, ++subOff) subIndices[subOff] = d;
    491:         }
    492:       }
    493:       PetscCall(ISCreateGeneral(comm, subSize, subIndices, PETSC_OWN_POINTER, &(*is)[i]));
    494:       /* Must have same blocksize on all procs (some might have no points) */
    495:       {
    496:         PetscInt bs = -1, bsLocal[2], bsMinMax[2];

    498:         bsLocal[0] = bs < 0 ? PETSC_INT_MAX : bs;
    499:         bsLocal[1] = bs;
    500:         PetscCall(PetscGlobalMinMaxInt(comm, bsLocal, bsMinMax));
    501:         if (bsMinMax[0] != bsMinMax[1]) {
    502:           bs = 1;
    503:         } else {
    504:           bs = bsMinMax[0];
    505:         }
    506:         PetscCall(ISSetBlockSize((*is)[i], bs));
    507:       }
    508:     }
    509:   }
    510:   /* Preserve discretizations */
    511:   if (len && dms[0]->probs) {
    512:     PetscCall(DMSetNumFields(*superdm, Nf));
    513:     for (i = 0, supf = 0; i < len; ++i) {
    514:       for (f = 0; f < Nfs[i]; ++f, ++supf) {
    515:         PetscObject disc;

    517:         PetscCall(DMGetField(dms[i], f, NULL, &disc));
    518:         PetscCall(DMSetField(*superdm, supf, NULL, disc));
    519:       }
    520:     }
    521:     PetscCall(DMCreateDS(*superdm));
    522:   }
    523:   /* Preserve nullspaces */
    524:   for (i = 0, supf = 0; i < len; ++i) {
    525:     for (f = 0; f < Nfs[i]; ++f, ++supf) {
    526:       (*superdm)->nullspaceConstructors[supf] = dms[i]->nullspaceConstructors[f];
    527:       if ((*superdm)->nullspaceConstructors[supf]) {
    528:         haveNull = PETSC_TRUE;
    529:         nullf    = supf;
    530:         oldf     = f;
    531:       }
    532:     }
    533:   }
    534:   /* Attach nullspace to IS */
    535:   if (haveNull && is) {
    536:     MatNullSpace nullSpace;

    538:     PetscCall((*(*superdm)->nullspaceConstructors[nullf])(*superdm, oldf, nullf, &nullSpace));
    539:     PetscCall(PetscObjectCompose((PetscObject)(*is)[nullf], "nullspace", (PetscObject)nullSpace));
    540:     PetscCall(MatNullSpaceDestroy(&nullSpace));
    541:   }
    542:   PetscCall(PetscSectionDestroy(&supersection));
    543:   PetscCall(PetscFree3(Nfs, sections, sectionGlobals));
    544:   PetscFunctionReturn(PETSC_SUCCESS);
    545: }