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 Note: 13 This should greatly improve the performance of the closure operations, at the cost of additional memory. 14 15 Level: intermediate 16 17 .seealso `DMPlexVecGetClosure()`, `DMPlexVecRestoreClosure()`, `DMPlexVecSetClosure()`, `DMPlexMatSetClosure()` 18 @*/ 19 PetscErrorCode DMPlexCreateClosureIndex(DM dm, PetscSection section) { 20 PetscSection closureSection; 21 IS closureIS; 22 PetscInt *clPoints; 23 PetscInt pStart, pEnd, sStart, sEnd, point, clSize; 24 25 PetscFunctionBegin; 26 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 27 if (!section) PetscCall(DMGetLocalSection(dm, §ion)); 28 PetscValidHeaderSpecific(section, PETSC_SECTION_CLASSID, 2); 29 PetscCall(PetscSectionGetChart(section, &sStart, &sEnd)); 30 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd)); 31 PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)section), &closureSection)); 32 PetscCall(PetscSectionSetChart(closureSection, pStart, pEnd)); 33 for (point = pStart; point < pEnd; ++point) { 34 PetscInt *points = NULL, numPoints, p, dof, cldof = 0; 35 36 PetscCall(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 PetscCall(PetscSectionGetDof(section, points[p], &dof)); 40 if (dof) cldof += 2; 41 } 42 } 43 PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points)); 44 PetscCall(PetscSectionSetDof(closureSection, point, cldof)); 45 } 46 PetscCall(PetscSectionSetUp(closureSection)); 47 PetscCall(PetscSectionGetStorageSize(closureSection, &clSize)); 48 PetscCall(PetscMalloc1(clSize, &clPoints)); 49 for (point = pStart; point < pEnd; ++point) { 50 PetscInt *points = NULL, numPoints, p, q, dof, cldof, cloff; 51 52 PetscCall(PetscSectionGetDof(closureSection, point, &cldof)); 53 PetscCall(PetscSectionGetOffset(closureSection, point, &cloff)); 54 PetscCall(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 PetscCall(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 PetscCall(DMPlexRestoreTransitiveClosure(dm, point, PETSC_TRUE, &numPoints, &points)); 66 PetscCheck(q * 2 == cldof, PetscObjectComm((PetscObject)dm), PETSC_ERR_PLIB, "Invalid size for closure %" PetscInt_FMT " should be %" PetscInt_FMT, q * 2, cldof); 67 } 68 PetscCall(ISCreateGeneral(PETSC_COMM_SELF, clSize, clPoints, PETSC_OWN_POINTER, &closureIS)); 69 PetscCall(PetscSectionSetClosureIndex(section, (PetscObject)dm, closureSection, closureIS)); 70 PetscCall(PetscSectionDestroy(&closureSection)); 71 PetscCall(ISDestroy(&closureIS)); 72 PetscFunctionReturn(0); 73 } 74