Actual source code: plexindices.c

  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: [](chapter_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;

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

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

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