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

  3: /*@
  4:   DMPlexCreateClosureIndex - Calculate an index for the given `PetscSection` for the closure operation on the `DM`

  6:   Not Collective

  8:   Input Parameters:
  9: + dm      - The `DM`
 10: - section - The section describing the layout in the local vector, or NULL to use the default section

 12:   Level: intermediate

 14:   Note:
 15:   This should greatly improve the performance of the closure operations, at the cost of additional memory.

 17: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSection`, `DMPlexVecGetClosure()`, `DMPlexVecRestoreClosure()`, `DMPlexVecSetClosure()`, `DMPlexMatSetClosure()`
 18: @*/
 19: PetscErrorCode DMPlexCreateClosureIndex(DM dm, PetscSection section)
 20: {
 21:   PetscSection closureSection;
 22:   IS           closureIS;
 23:   PetscInt    *clPoints;
 24:   PetscInt     pStart, pEnd, sStart, sEnd, point, clSize;

 26:   PetscFunctionBegin;
 28:   if (!section) PetscCall(DMGetLocalSection(dm, &section));
 30:   PetscCall(PetscSectionGetChart(section, &sStart, &sEnd));
 31:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
 32:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), &closureSection));
 33:   PetscCall(PetscSectionSetChart(closureSection, pStart, pEnd));
 34:   for (point = pStart; point < pEnd; ++point) {
 35:     PetscInt *points = NULL, numPoints, p, dof, cldof = 0;

 37:     PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points));
 38:     for (p = 0; p < numPoints * 2; p += 2) {
 39:       if ((points[p] >= sStart) && (points[p] < sEnd)) {
 40:         PetscCall(PetscSectionGetDof(section, points[p], &dof));
 41:         if (dof) cldof += 2;
 42:       }
 43:     }
 44:     PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points));
 45:     PetscCall(PetscSectionSetDof(closureSection, point, cldof));
 46:   }
 47:   PetscCall(PetscSectionSetUp(closureSection));
 48:   PetscCall(PetscSectionGetStorageSize(closureSection, &clSize));
 49:   PetscCall(PetscMalloc1(clSize, &clPoints));
 50:   for (point = pStart; point < pEnd; ++point) {
 51:     PetscInt *points = NULL, numPoints, p, q, dof, cldof, cloff;

 53:     PetscCall(PetscSectionGetDof(closureSection, point, &cldof));
 54:     PetscCall(PetscSectionGetOffset(closureSection, point, &cloff));
 55:     PetscCall(DMPlexGetTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points));
 56:     for (p = 0, q = 0; p < numPoints * 2; p += 2) {
 57:       if ((points[p] >= sStart) && (points[p] < sEnd)) {
 58:         PetscCall(PetscSectionGetDof(section, points[p], &dof));
 59:         if (dof) {
 60:           clPoints[cloff + q * 2]     = points[p];
 61:           clPoints[cloff + q * 2 + 1] = points[p + 1];
 62:           ++q;
 63:         }
 64:       }
 65:     }
 66:     PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points));
 67:     PetscCheck(q * 2 == cldof, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid size for closure %" PetscInt_FMT " should be %" PetscInt_FMT, q * 2, cldof);
 68:   }
 69:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPoints, PETSC_OWN_POINTER, &closureIS));
 70:   PetscCall(PetscSectionSetClosureIndex(section, (PetscObject)dm, closureSection, closureIS));
 71:   PetscCall(PetscSectionDestroy(&closureSection));
 72:   PetscCall(ISDestroy(&closureIS));
 73:   PetscFunctionReturn(PETSC_SUCCESS);
 74: }