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:   Note:
 13:   This should greatly improve the performance of the closure operations, at the cost of additional memory.

 15:   Level: intermediate

 17: .seealso 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;

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

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

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