1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 /*@ 4 DMPlexCreateClosureIndex - Calculate an index for the given `PetscSection` for the closure operation on the `DM` 5 6 Not Collective 7 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 11 12 Level: intermediate 13 14 Note: 15 This should greatly improve the performance of the closure operations, at the cost of additional memory. 16 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; 25 26 PetscFunctionBegin; 27 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 28 if (!section) PetscCall(DMGetLocalSection(dm, §ion)); 29 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 2); 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; 36 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; 52 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 } 75